{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homogeneous Gas Reactor (Ignition Zero D problem)\n",
    "\n",
    "***Energy equation***\n",
    "$$\n",
    "\\frac{dT}{dt}= -\\frac{1}{\\rho c_p}\\sum_{k=1}^{N_{spec}}\\dot{\\omega_{k}} W_k h_k = S_T\n",
    "$$\n",
    "***Species equation***\n",
    "$$\n",
    "\\frac{dY_k}{dt}=\\frac{1}{\\rho}\\dot{\\omega_{k}}W_k=S_{Y_k},\\,\\,\\,k=1\\ldots N_{spec}\n",
    "$$\n",
    "\n",
    "## Environment Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "TChem_install_directory ='where/TChem/is/installed'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append(TChem_install_directory+'/lib')\n",
    "import pytchem\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append(TChem_install_directory+'/example/runs/scripts/')\n",
    "import pmixSample\n",
    "import helper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pre-simulation\n",
    "* Set variables; temperature and stoichiometric ratio (fuel/air).\n",
    "* Convert from stoichiometric ratio to mass fraction of H2 and O2.\n",
    "* Create samples spaning over the variable ranges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pressure, Temperature, phi(stoichiometric ratio)\n",
    "one_atm = 101325 # [Pa]\n",
    "Temp = 900   #  temperature [K]\n",
    "Pressure = 1*one_atm; # [Pa]\n",
    "phi = 1.0; # Maximum equivalence ratio [-]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=100\n",
    "Nvar = 6\n",
    "sample = np.zeros([N,Nvar])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fuel =\"H2\"\n",
    "header_fuel = \"T P \" + fuel + \" O2\"\n",
    "Nvar = 4\n",
    "Yp_fuel, Yr_o2 = pmixSample.getMassFractionH2(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = np.zeros([N,Nvar])\n",
    "sample[:,0] = Temp\n",
    "sample[:,1] = Pressure\n",
    "sample[:,2] = Yp_fuel\n",
    "sample[:,3] = Yr_o2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TChem Simulation\n",
    "\n",
    "### Initialize TChem Driver Object\n",
    "\n",
    "* Initialization of Kokkos.\n",
    "* Create a TChem driver object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pytchem.initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "tchem = pytchem.TChemDriver()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Get help from TChem driver object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on TChemDriver in module pytchem object:\n",
      "\n",
      "class TChemDriver(pybind11_builtins.pybind11_object)\n",
      " |  A class to manage data movement between numpy to kokkos views in TChem::Driver object\n",
      " |  \n",
      " |  Method resolution order:\n",
      " |      TChemDriver\n",
      " |      pybind11_builtins.pybind11_object\n",
      " |      builtins.object\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __init__(...)\n",
      " |      __init__(self: pytchem.TChemDriver) -> None\n",
      " |  \n",
      " |  cloneGasKineticModel(...)\n",
      " |      cloneGasKineticModel(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Internally create clones of the kinetic model\n",
      " |  \n",
      " |  computeGasEnthapyMass(...)\n",
      " |      computeGasEnthapyMass(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute enthalpy mass and mixture enthalpy\n",
      " |  \n",
      " |  computeGasNetProductionRatePerMass(...)\n",
      " |      computeGasNetProductionRatePerMass(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute net production rate\n",
      " |  \n",
      " |  computeGasReactionRateConstants(...)\n",
      " |      computeGasReactionRateConstants(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute forward/reverse rate constant\n",
      " |  \n",
      " |  computeJacobianHomogeneousGasReactor(...)\n",
      " |      computeJacobianHomogeneousGasReactor(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute Jacobian matrix for homogeneous gas reactor\n",
      " |  \n",
      " |  computeRHS_HomogeneousGasReactor(...)\n",
      " |      computeRHS_HomogeneousGasReactor(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Compute RHS for homogeneous gas reactor\n",
      " |  \n",
      " |  computeTimeAdvanceHomogeneousGasReactor(...)\n",
      " |      computeTimeAdvanceHomogeneousGasReactor(self: pytchem.TChemDriver) -> float\n",
      " |      \n",
      " |      Compute Time Advance for a Homogeneous-Gas Reactor\n",
      " |  \n",
      " |  createAllViews(...)\n",
      " |      createAllViews(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate all necessary workspace for this driver\n",
      " |  \n",
      " |  createGasEnthapyMass(...)\n",
      " |      createGasEnthapyMass(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for enthalpy mass  (# samples, # species )\n",
      " |  \n",
      " |  createGasKineticModel(...)\n",
      " |      createGasKineticModel(self: pytchem.TChemDriver, chemkin_input: str, thermo_data: str) -> None\n",
      " |      \n",
      " |      Create a kinetic model from CHEMKIN input files\n",
      " |  \n",
      " |  createGasKineticModelConstData(...)\n",
      " |      createGasKineticModelConstData(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Internally construct const object of the kinetic model and load them to device\n",
      " |  \n",
      " |  createGasKineticModelConstDataWithArreniusForwardParameters(...)\n",
      " |      createGasKineticModelConstDataWithArreniusForwardParameters(self: pytchem.TChemDriver, reac_indices: numpy.ndarray[numpy.int32], factors: numpy.ndarray[numpy.float64]) -> None\n",
      " |      \n",
      " |      Creates clones of the kinetic model; modifies the Arrhenius forward parameters of the clones;and creates a const object of the kinetics models.factors is a 3d array of size: number of samples, number of reactions to be modified, and 3 kinetic parameters. kinetic parameters: pre exponential (0), temperature coefficient(1), and activation energy(2)\n",
      " |  \n",
      " |  createGasNetProductionRatePerMass(...)\n",
      " |      createGasNetProductionRatePerMass(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for net production rate per mass (# samples, # species)\n",
      " |  \n",
      " |  createGasReactionRateConstants(...)\n",
      " |      createGasReactionRateConstants(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for forward/reverse rate constants  (# samples, # reactions )\n",
      " |  \n",
      " |  createJacobianHomogeneousGasReactor(...)\n",
      " |      createJacobianHomogeneousGasReactor(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for homogeneous-gas-reactor Jacobian  (# samples, # species + 1  # species + 1)\n",
      " |  \n",
      " |  createRHS_HomogeneousGasReactor(...)\n",
      " |      createRHS_HomogeneousGasReactor(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for homogeneous-gas-reactor RHS  (# samples, # species + 1 )\n",
      " |  \n",
      " |  createStateVector(...)\n",
      " |      createStateVector(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Allocate memory for state vector (# samples, state vector length)\n",
      " |  \n",
      " |  freeAllViews(...)\n",
      " |      freeAllViews(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Free all necessary workspace for this driver\n",
      " |  \n",
      " |  getGasArrheniusForwardParameter(...)\n",
      " |      getGasArrheniusForwardParameter(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getGasArrheniusForwardParameter(self: pytchem.TChemDriver, reac_indices: numpy.ndarray[numpy.int32], param_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive pre exponential for reactions listed by reaction_indices\n",
      " |      \n",
      " |      2. getGasArrheniusForwardParameter(self: pytchem.TChemDriver, imodel: int, reac_indices: numpy.ndarray[numpy.int32], param_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive pre exponential for reactions listed by reaction_indices\n",
      " |  \n",
      " |  getGasEnthapyMass(...)\n",
      " |      getGasEnthapyMass(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive enthalpy mass per species for all samples\n",
      " |  \n",
      " |  getGasEnthapyMixMass(...)\n",
      " |      getGasEnthapyMixMass(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve mixture enthalpy for all samples\n",
      " |  \n",
      " |  getGasForwardReactionRateConstants(...)\n",
      " |      getGasForwardReactionRateConstants(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getGasForwardReactionRateConstants(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive forward rate constants for a single sample\n",
      " |      \n",
      " |      2. getGasForwardReactionRateConstants(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve forward rate constants  for all samples\n",
      " |  \n",
      " |  getGasNetProductionRatePerMass(...)\n",
      " |      getGasNetProductionRatePerMass(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getGasNetProductionRatePerMass(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive net production rate for a single sample\n",
      " |      \n",
      " |      2. getGasNetProductionRatePerMass(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve net production rate for all samples\n",
      " |  \n",
      " |  getGasReverseReactionRateConstants(...)\n",
      " |      getGasReverseReactionRateConstants(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getGasReverseReactionRateConstants(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive reverse rate constants for a single sample\n",
      " |      \n",
      " |      2. getGasReverseReactionRateConstants(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve reverse rate constants  for all samples\n",
      " |  \n",
      " |  getJacobianHomogeneousGasReactor(...)\n",
      " |      getJacobianHomogeneousGasReactor(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive homogeneous-gas-reactor Jacobian for a single sample\n",
      " |  \n",
      " |  getLengthOfStateVector(...)\n",
      " |      getLengthOfStateVector(self: pytchem.TChemDriver) -> int\n",
      " |      \n",
      " |      Get the size of state vector i.e., rho, P, T, Y_{0-Nspec-1}\n",
      " |  \n",
      " |  getNumberOfReactions(...)\n",
      " |      getNumberOfReactions(self: pytchem.TChemDriver) -> int\n",
      " |      \n",
      " |      Get the number of reactions registered in the kinetic model\n",
      " |  \n",
      " |  getNumberOfSamples(...)\n",
      " |      getNumberOfSamples(self: pytchem.TChemDriver) -> int\n",
      " |      \n",
      " |      Get the number of samples which is currently used in the driver\n",
      " |  \n",
      " |  getNumberOfSpecies(...)\n",
      " |      getNumberOfSpecies(self: pytchem.TChemDriver) -> int\n",
      " |      \n",
      " |      Get the number of species registered in the kinetic model\n",
      " |  \n",
      " |  getRHS_HomogeneousGasReactor(...)\n",
      " |      getRHS_HomogeneousGasReactor(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getRHS_HomogeneousGasReactor(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrive homogeneous-gas-reactor RHS for a single sample\n",
      " |      \n",
      " |      2. getRHS_HomogeneousGasReactor(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve homogeneous-gas-reactor RHS_ for all samples\n",
      " |  \n",
      " |  getSpeciesIndex(...)\n",
      " |      getSpeciesIndex(self: pytchem.TChemDriver, species_name: str) -> int\n",
      " |      \n",
      " |      Get species index\n",
      " |  \n",
      " |  getStateVariableIndex(...)\n",
      " |      getStateVariableIndex(self: pytchem.TChemDriver, var_name: str) -> int\n",
      " |      \n",
      " |      Get state variable index\n",
      " |  \n",
      " |  getStateVector(...)\n",
      " |      getStateVector(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. getStateVector(self: pytchem.TChemDriver, sample_index: int) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve state vector for a single sample\n",
      " |      \n",
      " |      2. getStateVector(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve state vector for all samples\n",
      " |  \n",
      " |  getTimeStep(...)\n",
      " |      getTimeStep(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve time line of all samples\n",
      " |  \n",
      " |  getTimeStepSize(...)\n",
      " |      getTimeStepSize(self: pytchem.TChemDriver) -> numpy.ndarray[numpy.float64]\n",
      " |      \n",
      " |      Retrieve time step sizes of all samples\n",
      " |  \n",
      " |  modifyGasArrheniusForwardParameters(...)\n",
      " |      modifyGasArrheniusForwardParameters(self: pytchem.TChemDriver, reac_indices: numpy.ndarray[numpy.int32], factors: numpy.ndarray[numpy.float64]) -> None\n",
      " |      \n",
      " |      Modify the cloned kinetic models Arrhenius parameters\n",
      " |  \n",
      " |  setNumberOfSamples(...)\n",
      " |      setNumberOfSamples(self: pytchem.TChemDriver, number_of_samples: int) -> None\n",
      " |      \n",
      " |      Set the number of samples; this is used for Kokkos view allocation\n",
      " |  \n",
      " |  setStateVector(...)\n",
      " |      setStateVector(*args, **kwargs)\n",
      " |      Overloaded function.\n",
      " |      \n",
      " |      1. setStateVector(self: pytchem.TChemDriver, sample_index: int, 1d_state_vector: numpy.ndarray[numpy.float64]) -> None\n",
      " |      \n",
      " |      Overwrite state vector for a single sample\n",
      " |      \n",
      " |      2. setStateVector(self: pytchem.TChemDriver, 2d_state_vector: numpy.ndarray[numpy.float64]) -> None\n",
      " |      \n",
      " |      Overwrite state vector for all samples\n",
      " |  \n",
      " |  setTimeAdvanceHomogeneousGasReactor(...)\n",
      " |      setTimeAdvanceHomogeneousGasReactor(self: pytchem.TChemDriver, tbeg: float, tend: float, dtmin: float, dtmax: float, jacobian_interval: int, max_num_newton_iterations: int, num_time_iterations_per_interval: int, atol_newton: float, rtol_newton: float, atol_time: float, rtol_time: float) -> None\n",
      " |      \n",
      " |      Set time advance object for homogeneous gas reactor\n",
      " |  \n",
      " |  showViewStatus(...)\n",
      " |      showViewStatus(self: pytchem.TChemDriver) -> None\n",
      " |      \n",
      " |      Print member variable view status\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Static methods inherited from pybind11_builtins.pybind11_object:\n",
      " |  \n",
      " |  __new__(*args, **kwargs) from pybind11_builtins.pybind11_type\n",
      " |      Create and return a new object.  See help(type) for accurate signature.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(tchem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create Kinetic Model \n",
    "\n",
    "* Inputs are the reactions mechanism files; in this case, we use the GRI3.0 gas reaction mechanism"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs_directory = TChem_install_directory + '/example/data/H2/'\n",
    "tchem.createGasKineticModel(inputs_directory+'chem.inp',inputs_directory+'therm.dat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "tchem.setNumberOfSamples(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### Model Variation\n",
    "* Set list of reaction to be modified ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18]\n"
     ]
    }
   ],
   "source": [
    "Nparameter= 19\n",
    "reaction_index = np.arange(Nparameter)\n",
    "print(reaction_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Set a 3d array, factors, of dimension: Number of samples (N), number of reaction to be modified (Nparameter), 3 that corresponds to three Arrenius parameters: pre exponetial (0), temperature coeficient(1), activation energy(2). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "μ=0 \n",
    "σ=.25\n",
    "A_modifier = np.random.lognormal(μ, σ, size=(N, Nparameter))\n",
    "factors = np.empty([N,Nparameter,3])\n",
    "ones = np.ones([N,Nparameter])\n",
    "factors[:,:,0] = A_modifier # //pre exponetial \n",
    "factors[:,:,1] = ones # temperature coef \n",
    "factors[:,:,2] = ones # activation energy "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Create a list of const objects of the kinetic model with the modified parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "tchem.createGasKineticModelConstDataWithArreniusForwardParameters(reaction_index,factors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Verify the changes comparing reference values vs modified variables for a specific sample (isample)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "ver_reaction_index = np.zeros(3,dtype=int)\n",
    "ver_reaction_index[0] = 1\n",
    "ver_reaction_index[1] = 3\n",
    "ver_reaction_index[2] = 6\n",
    "pre_exponetial_index = 0\n",
    "isample=84\n",
    "\n",
    "A_original = tchem.getGasArrheniusForwardParameter(ver_reaction_index, pre_exponetial_index); \n",
    "A_modified = tchem.getGasArrheniusForwardParameter(isample, ver_reaction_index, pre_exponetial_index);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pre exponential Original: [5.080e+04 1.230e+04 4.714e+18]\n",
      "Pre exponential Modified: [6.34132474e+04 1.51141650e+04 5.65527868e+18]\n",
      "Factor                  : [1.24829227 1.2287939  1.19967728]\n"
     ]
    }
   ],
   "source": [
    "print('Pre exponential Original:', A_original)\n",
    "print('Pre exponential Modified:', A_modified)\n",
    "print('Factor                  :', A_modified/A_original)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "for isample in [3, 20, 40, 53, 83]:\n",
    "    A_original = tchem.getGasArrheniusForwardParameter(reaction_index, pre_exponetial_index); \n",
    "    A_modified = tchem.getGasArrheniusForwardParameter(isample, reaction_index, pre_exponetial_index);\n",
    "    factor = A_modified/A_original\n",
    "    for i, ifac in enumerate(factor):\n",
    "        error = (ifac - A_modifier[isample,i])/ A_modifier[isample,i]\n",
    "        if abs(error) > 1e-15:\n",
    "            print('error in ',error,' in sample No',isample,' factor No', i)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set State Vector\n",
    "\n",
    "* Get index for variables. \n",
    "* Pass a numpy array to the TChem object to set the state vector. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "Variables = ['T','P','H2','O2']\n",
    "indx=[]\n",
    "for var in Variables:\n",
    "    indx += [tchem.getStateVariableIndex(var)]\n",
    "\n",
    "state = np.zeros([N, tchem.getLengthOfStateVector()])\n",
    "for sp in range(N):\n",
    "    state[sp,indx] = sample[sp,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "tchem.createStateVector()\n",
    "tchem.setStateVector(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Time Integrator and Its Parameters \n",
    "* Kokkos team policy is constructed for parallel execution\n",
    "* Workspace per team is also allocated with the teampolicy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "tend = 2\n",
    "tchem.setTimeAdvanceHomogeneousGasReactor(tbeg=0, \n",
    "                     tend=tend, \n",
    "                     dtmin=1e-10,\n",
    "                     dtmax=1e-3, \n",
    "                     atol_time=1e-12,\n",
    "                     rtol_time=1e-6,\n",
    "                     max_num_newton_iterations=20,\n",
    "                     atol_newton=1e-18,\n",
    "                     rtol_newton=1e-8,\n",
    "                     num_time_iterations_per_interval=20,\n",
    "                     jacobian_interval=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Advance simulations: for each iteration, get state vector. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count  50 simulation time  4.450e-01 s, total time  2.000e+00 s\n",
      "count 100 simulation time  1.445e+00 s, total time  2.000e+00 s\n",
      "Wall time is 0.3973 mins for time integration\n"
     ]
    }
   ],
   "source": [
    "solution = []\n",
    "time_simulation = []\n",
    "\n",
    "tavg = 0\n",
    "t0 = time.time()\n",
    "icount = 0\n",
    "while (tavg < tend*0.999):\n",
    "    tavg = tchem.computeTimeAdvanceHomogeneousGasReactor()\n",
    "    solution += [tchem.getStateVector()]\n",
    "    time_simulation += [tchem.getTimeStep()]\n",
    "    icount +=1\n",
    "    if (icount%50 == 0):\n",
    "        print(f\"count {icount:3d} simulation time {tavg:10.3e} s, total time {tend:10.3e} s\")\n",
    "\n",
    "t1 = time.time()\n",
    "\n",
    "print(f\"Wall time is {(t1 - t0)/60:0.4f} mins for time integration\")\n",
    "\n",
    "solution = np.array(solution)\n",
    "time_simulation = np.array(time_simulation) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post-processing\n",
    "* Get variable index for plotting. \n",
    "* Compute the ignition delay time with data produced by the simulation.\n",
    "* Plot time profiles for important variables and for the ignition delay time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_indx   = tchem.getStateVariableIndex('T')\n",
    "H2_indx = tchem.getStateVariableIndex('H2')\n",
    "O2_indx  = tchem.getStateVariableIndex('O2')\n",
    "OH_indx  = tchem.getStateVariableIndex('OH')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_threshold=1500\n",
    "def computeIgnitionDelayTime(T_threshold=1500):\n",
    "    ntimes, nsamples, nVars = np.shape(solution)\n",
    "    ignDelayTime = []\n",
    "    for isample in range(nsamples):\n",
    "        Indx = np.where(solution[:,isample,T_indx] >= T_threshold )\n",
    "        ignDelayTime += [time_simulation[Indx[0][0],isample] ] \n",
    "    return np.array(ignDelayTime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "IgnTime = computeIgnitionDelayTime()\n",
    "np.savetxt('QoI_H2.txt',IgnTime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAEcCAYAAACYg/MAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABFwUlEQVR4nO3deXxU5dn/8c81CWGLGwIFFQSUx6WiqEAbazGKFteftra17rgraMWtfdxxqWul2Kp9irZFaam0tYtiRVmMWo2lgKggWESgiiIiIrKGJNfvj/skTIZJMgmzJXzfvs4rZ7nnnGvGE3LNfe7F3B0RERERyU+xXAcgIiIiIvVTsiYiIiKSx5SsiYiIiOQxJWsiIiIieUzJmoiIiEgeK8x1APkoFot5+/btcx2GiEirtH79end3VRaIpEjJWhLt27dn3bp1uQ5DRKRVMrMNuY5BpCXRNxsRERGRPKZkTURERCSPKVkTERERyWNK1kRERETymJI1ERERkTymZE1EREQkjylZExEREcljStZERERE8pgGxU1iww4b6PNgH2IWw8yIWSyss2W9bWFbiouKKS4qZoeiHWp/9typJ/t32Z/9u+zPhhUbOOmkk1i/fn2u35KIiIi0UErWkohVxTi85+FUezWOh5/utdtV1VVsqtrE2oq1LFuzjLUVa/my4kvWbFrD+s1bErN9O+7LgnULOPGgE+ncuXMO35GISP4YN25crkMQaVHM3XMdQ97p2LGjN2e6KXdnxboVzF85n38v+zf3vXQfK1ev5JVTX+Hwrx2egUhFRFoeM1vv7h1zHYdIS6GatTQyM75S/BW+UvwVSnuV8sl7n/DA5gf4aMNHuQ5NREREWih1MMigNtYGgCqvynEkIiIi0lIpWcsgM8t1CCIiItLCKVkTERERyWNK1kRERETymJI1ERERkTymZE1EREQkjylZExEREcljStZERERE8piSNREREZE8pmRNREREJI8pWRMRERHJY0rWRERERPKYkjURERGRPKZkTURERCSPKVkTEZG8ZGbDzWyxmW00s1lm9s1Gyg81s3Iz+9LMVprZ383sf7IVr0imFOY6ABERkURmdhrwIDAc+Gf08zkz29/d/5ukfG/g78DPgbOBYuA+4B/A3tmKO1WxWMzbt2+f6zAkx9avX+/u3mjFmZI1ERHJR1cD49z90Wj7CjM7FrgMuD5J+UOBNsD17l4FYGZ3A9PNrLO7r8xG0Klq374969aty3UYkmNmtiGVcnoMKiIiecXMigjJ1wsJh14ADqvnZTOBzcCFZlZgZjsA5wL/zrdETaSplKyJiEi2FZrZzLjl4oTjnYEC4JOE/Z8A3ZKd0N2XAMcAtwGbgC+AfsCJ6Qxc0qe8vJy7776b8vLyXIeS9/QYVEREsq3S3QekUM4Tti3JvnDArBvwa+AJ4A/ADsDtwB/N7Ch3r96GeFuV8g/KKVtSRmmvUkp6lGTtOvHbbz/3NiNGjKCqqoo2bdpw2223se+++2YslpZOyZqIiOSblUAVW9eidWXr2rYaI4B17v6jmh1mdhbwAeHR6T8zEGezbey4kaOfOBozw7CtfgJpP2ZmfLruU6a+P5Uqr6LACjih7wnstsNuFMQKKLCCrX4WxgrrPVYQi44nObbo80Xc8fIdVFZVUlhQyCWHXMLY2WPZXLUZx/EKhx8CZVAxu4Lrr0/WDFFqKFkTEZG84u4VZjaL8FjzT3GHjgGequdlHQgJXrya7fxr8mOwoXID7h6Sl7ifwFb7tvWYRxWSn2/4nKrQ/4Iqr+LFJS/SrrAdVV5FVXVV0p/bqqKqgl/8+xd1d7aJlpMgFovxm8t/w0EHHbTN12ppDj744JTKKVkTEZF8NBoYb2YzgFeBS4HdgP+D2p6eg9x9SFT+WeAqM7sVmEB4DHoXoWZtVpZjb1S7te149fxXs37d8g/KGfLEECqqKigqKOL5s55v8FGou1Pt1fUmc5XVlUmPzfpoFhc8fQGV1ZUUxgoZ+fWRjHl9DJuqNoUTW80FoOcpPTn33HMz/+ZbMCVrIiKSd9x9opntCtwEdAfmAse7+9KoSHdgr7jy083sDOBHwHXABuB14Fh31xgZkZIeJUw7Z1rKbdbMLDzepCB0+UjRvp33pc8ufepc5+R9Tmb4s8OZ88mc0PLQwtJz957b8pa2C0rWREQkL7n7I8Aj9RwblmTfk8CTGQ6rxSvpUZLRjgX1XaekRwlvXPoGe4zeg2VrltXu79SuU8Zjaeny7zm+iIiItFoDdxq4ZcMB1Xs2SsmaiIiIZM/H0U9P2JZ6KVkTERGRrOnWOxqRJRo1r3Zb6qVkTURERLLmnMHnhBWHNgVttmxLvZSsiYiISE7EYkpDUqFPSURERLKmbElZWDGorK7csi31UrImIiIiWVPaqzSsOBQVFG3ZlnopWRMREZGsKelRQtHmIjp80YFp50zLyphvLZ2SNREREcmqmMdo/2V7JWopUrImIiIiWVVt1WzYcQPlH5TnOpQWQcmaiIiIZE35B+VUtKlg/Y7rGfLEECVsKVCyJiIiIlkT3xu0oqpCvUFToGRNREREsia+N2ihFao3aAqUrImIiEj2fAh8CSwHf9zDtjRIyZqIiIhkTVlZGVQCn0LVkqqwLQ0qrO+Amf23Gedz4AR3n9v8kERERKS1Ki0thSfDelFRUdiWBtWbrAF7AP8APk3xXDHgLKBoW4MSERGR1qmkpIS2f21LQYcCpk6bSkmJxlprTEPJGsDt7j4jlROZWSFw9raH1Hq4e65DEBERyTuxWIwOHTooUUtRQ23WbgY+SPVE7l4ZvWbZtgYlIiIirVdVrIr1xes1xlqK6k3W3P0n7v5xKicxs+K413ySruBaC8NyHYKIiEheKP+gnIqiCtbvoEFxU1VvsmZmF6RyAjPbGZiSroBaI0ePQ0VERECD4jZHQ49Bx5rZmQ292Mw6AdOAA9MaVSuhGjUREZG64gfFLSoo0qC4KWiog8HvgN+a2UZ3fyrxoJl1ISRqfYCTMxSfiIhIxplZDBgE9ATaJR539yeyHlQrVdKjhKJNRbTZ3IYpV06hpIc6GTSmoWRtGGEYjglm9l13f6bmgJl1JyRqewDHu/vLGY1SREQkQ8xsf+BvwF6Q9JGIA0rW0qiguoD269orUUtRvcmau7uZnQX8EfijmZ3s7i+YWU9CotYFGOruahkoIiIt2SOEv4ffB94GNuU2nO2DmgqlrsHppty9CjgNmAr81cyGAS8DnYAhStRERKQVOAS41t2fcvf/uPvSxCXXAbY21dXVrF+/nvJypRGpaHRu0Gj8tFOBV4BfA+2BUnefleHYREREsmElUJHrILYX5eXlbNq4iXXr1zFkyBAlbCloaG7QxOfzXxKe5b8PXGdWp/rS3f3c9IcnIiKScT8DRpjZc9ETJcmg+InbKyoqKCsr00wGjWiog8Fg2GqAsKVAt2iJp4HERESkpeoC7AO8Y2ZTgFUJx93db81+WK1TaWkpTAjrmsg9NQ11MOiV7ouZ2fXAdwi/FJuA14Hr3X1uXJlxQGIt3b/c/etxZdoCPwVOJzyWnQYMd/cP48rsAvwc+H/RrqeBK9x9dXrflYiItHA3xa33TXLcASVraVJSUkLbp9pS2LGQKdOmqFYtBY22WUuzUkKvm8OAo4BKYGo0uG68qUD3uOX4hONjCO3oTge+CewITDKzgrgyEwiNRo8Djo3Wx6fvrYiISGvg7rFGloLGzyJNEYvF6NBRE7mnqqHHoGnn7kPjt83sbOAL4BvAM3GHNrn78mTnMLOdgAuA89x9Stx5lgJHA8+b2X6EBO1wd38tKnMJ8IqZ7ePu76b3nYmIiEhTaOiO1DU0N2iVmQ1K9URmVhC95pAmXH+HKIbPE/YfbmYrzOw/ZvaomXWNO3Yo0AZ4oWaHu38AzCfU2AGUAGuB1+Je9yqwLq5MYvwXm9lMM5tZWVnZhLcgIiKtgZmdaGb3m9mvzew+Mzsh1zG1Zpo3O3UNPQY1oJuZ9UxlAfYk+cjPDXkQmAPE99udDJwDDAGuIUz/MT1qpwahc0MVoat1vE/Y0vGhG/Cpu9feCdH6CrbuHFFzfKy7D3D3AYWFWa1wFGk1zKzBZdiwYbkOMe3GjRtHcXFxrsOQbWBmO5jZS4S2zVcSmt6MBJ42szIz0//gdPPw74WkprGs5K9NPF/KabKZjQYOJzyqrO0q7e5PxhV728xmER5xngD8paFTJlw/WSyJZUQkjT7++OPa9UmTJnHRRRfV2de+fftchNUsFRUVFBUVtfprCgB3Edo1nw086e5VURvoHwC/jI7/MIfxyXauoZq184Dzm7EsbuyiZvYzQueAo9z9/YbKuvtHwIds6aGzHCgAOicU7UqoXasp09Xi0vZovUtcGRFJs27dutUuO++881b7Xn75ZQ499FDatWtH7969ufHGG6mo2DIWaa9evbj99tsZNmwYO+ywAz169GDixImsXr2aH/zgBxQXF9O3b19eeKG2FQRlZWWYGZMmTaJ///60a9eOQw89lFmz6o7b/dprr3HEEUfQoUMHdt99dy677DLWrFlTe7y0tJTLLruMa6+9li5duvCNb3wDgNGjR3PggQfSsWNHdt99dy688EJWr15de+3zzjuPdevW1dYejho1qva9/PSnP60TQ2lpKZdffnmd9ztq1CjOP/98dt55Z84888yUYpW0OxW4yd1/X1N54O5V7v574ObouKRRVayKde3XUf6BBsRNRb3Jmrs/3swlsf1ZHWb2IHAGIVFb0FiAZtYZ2B2o+Xo+C9gMHBNXZg9gP7a0USsHiglt12qUAB2p245NRLLk+eef58wzz+Tyyy9n3rx5/OY3v+HPf/4zN9xwQ51yY8aMYdCgQcyePZvvf//7nHvuuZxxxhkcf/zxzJkzh8GDB3PWWWexcePGOq+79tpruffee5k5cyZ9+vThhBNOYP369QC8/fbbfOtb3+L//b//x5tvvslf/vIX5syZw/nnn1/nHL/73e9wd1555RWeeCKMCx6LxRgzZgzz5s1jwoQJzJgxgyuuuAKAww47jDFjxtChQwc+/vhjPv74Y6699tomfS6jR49m3333ZebMmdx1110px7o9MLPhZrbYzDaa2Swz+2Yj5c3MRprZAjPbZGYfm9k9KVxqV+Cdeo69Ex2XNCn/oJyK9hWs67iOIU8MUcKWgqw2zjKzhwnVzKcAn5tZTfuxte6+NmoXMAp4ipCc9QLuJrQ1+yuAu39hZr8G7jezFcBnwGjgLcKQH7j7fDObDPzKzC4iPP78FTBJPUGlJRs5ciRz5szJ6jX79+/PmDFjtvk8P/nJT7juuus477zzANhrr7249957Oeuss7j//vtr268MHTqU4cOHA3DbbbcxevRo9t57b8455xwAbr75Zn7zm98wd+5cBgwYUHv+m2++maFDQ4fz3/72t+yxxx5MmDCBCy+8kPvvv5/TTjuNa665prb8L3/5Sw4++GBWrFhB166hD1Pv3r154IEH6sQ9cuTI2vVevXpx3333cfLJJ/P4449TVFTETjvthJnRrVvS5rCNOuKII/jRj35Uu33OOeekFGtrZ2anEdo1Dwf+Gf18zsz2d/f/1vOyB4ATgesIE7LvRBj+qTGLo9dNSXLseFJ4YiSpK1tSFlYMKqoqKFtSRkkPDeHRkGy3pB8e/ZyWsP82QpJWBfQjdDDYmZCwvQh8392/jCt/FWGMtolsGRT3nIRpQs4kDIpb87zkaeByRCQnZs2axYwZM7j33ntr91VXV7NhwwaWL19O9+7hb+qBBx5Ye7y4uJgOHTrQr1+/2n1f+cpXAFixYkWd88eP11RcXEy/fv145513aq/93nvvMXHixNoyNf2PFi1aVJsAHXrooVvFPX36dO6++27mz5/PF198QVVVFRUVFSxfvpzddtuteR9GnPiEsymxbgeuBsa5+6PR9hVmdixwGXB9YmEz2we4AjjQ3efHHXojhWv9CnggqjD4PeFvTzdCm7ULo1gkTUp7lYYVh6LCoi3bUq9sj7PWYNcPd98ADG2oTFRuI+GX8ooGyqwCzmpqjCL5LB01XLlSXV3Nrbfeyve+972tjnXp0qV2vU2bNnWOmVmdfTU1cNXV1U269oUXXshVV1211bHdd9+9dr1jx451ji1dupQTTjiBiy66iNtvv51dd92V2bNnc/rpp9dpa5dMLBarTbJqbN68eatyiddMNdbWzMyKCMM0/TTh0AvUM/wScDJh7upjzexZQjOfl4Dr3H1FPa8BwN1/ZmZdCBUBw2rCIMy0c4+7P9ic9yHJlfQooc26NrAZxgwdo1q1FGiMChHJikMOOYQFCxaw9957Z+T8r7/+On369AFg3bp1zJ07t/bR6SGHHMK8efOafO2ZM2dSUVHBz372MwoKwiD2kyZNqlOmqKiIqqqt5/7u0qVLnZ6wGzduZMGCBRx88MENXrO5sbYwhWY2M257rLuPjdvuTOhIltgh7BPC4OfJ9CEMIfUDQsLlhGTvGTMrcfcGs3t3v8HM7ge+DnQizA/6emPtsKXpysvL2bx+M6yEkaeOpN+0fprJoBHZnm5KRLZTt9xyCxMmTOCWW25h7ty5LFiwgD//+c912mttizvvvJMpU6Ywb948zj//fIqKijjjjDMA+PGPf8yMGTO49NJLeeONN3jvvfeYNGkSl1xySYPn7Nu3L9XV1YwZM4bFixfzhz/8YavazV69erFx40amTJnCypUrazs1HHXUUfz+97+nrKysNqZkNWuJmhtrC1NZM65ltIytp1ziUEsNDb8UA9oCZ7v7y+7+CqGN9CBgYCpBufvn7v5c1Cv0OSVqmVFWVla7XlFRUWdbkmtSsmZmMTM7wMyOMLOOjb9CRCQYOnQozz77LC+++CKDBg1i0KBB3HPPPfTs2TMt57/nnnu45pprOOSQQ1i4cCGTJk2qfcR44IEH8vLLL7NkyRKOOOIIDjroIK6//vra9m/1OfDAA3nwwQcZPXo0+++/P4899thWw3EcdthhXHrppZx++ul06dKF++67D4Drr7+eo446ipNPPplvfetbHH744RxySOMTvDQ31lZmJaENc2KvjfghmhJ9TEgC/xO3byGhffNWN5mZDa4Z7DZab3DZ5ncktUpLS2tT7qKiorAtDbLENhX1FjQbAdzKli7MA919tpn9DZju7j/PTIjZ17FjR1+3bt02n+eGJ27g7sV3M2HwBE4/8vQ0RCYiicrKyjjyyCP59NNP6dw5cfhFyUdmtt7dG/zCb2b/At5094vj9v0HeMrdk3Uw+BbwPLC3uy+K9u0FvAd8zd1nJJSvBr7u7jOi9fr+GBphEpy0Tuaerr8zLVW7q9pRuKaQKRdO2a4fgabyuwAptlmLhr94EPgNoYHnH+MOv0IYMLDVJGsiIpJzo4HxZjaDMLfzpcBuwP8BmNndwCB3HxKVnwrMBn5jZiOjfWOAfwHx7eNqHMmWsdWOQrPbZFWsIEZxcfF2nag1RaodDK4GHnD3H0dTcMRbQBjTRkREJC3cfaKZ7QrcRBgrbS5wvLsvjYp0B/aKK19tZicSKg5eBjYQxk27OlnnAnd/KW69LFPvQ+qh1LhJUk3WehOql5NZRxgTTUQk60pLS7caIkNaB3d/BHiknmPDkuz7GNh6bJhGmNn7wLfd/c0kxw4Annb3Pk09r0i6pNrBYCVhNoFk9gGWpSUaERGR7OtF6EmaTDvCkCAiOZNqsvYMcIuZxX+z8GjezquAv6U7MBERkSyqr3p2ALA6i3GIbCXVZO0mwkjOcwmNOJ3QLmA+oXv17RmJTkREJAPM7Coz+6+Z/ZfwN+2Zmu245VPgYWBybqNtfaqrqlm7bi3l5ZrEPRUptVlz98/MbAAwkjAd1KLotQ8BP3P3NRmLUEREJP3eZ8s81ecSeox+mlBmE6HH6GNZjKvVKy8vZ9OmTWxat4khQ4Ywbdo09QptRKPJWjRH273ABHe/A7gj41GJiIhkkLv/Hfg71M43e7u7L85pUNuJZDMYKFlrWKOPQd29ArgEaJ/5cFonVx9lEZF8dgmQdLJ3M+toZm2yHE+rFj9jgWYwSE2qbdbeAPplMhAREZEceTRakvlVtEialJSU0LZtW4qLi/UINEWpJmvXANea2YkW1RdL4wx9VLL9GTZsGCeeeOJW+2fOnImZsWTJEt58801OP/10evToQfv27dlnn324//77qa7eauxSkWw4kuiRaBJPA0PqOSbNFLMYHYs7KlFLUaqD4v4J2IlwM1ea2QrqdnN2d9c4NCKSklmzZtGlSxfGjx9Pz549mTFjBhdddBGbN2/mhhtuyHV4sv3pSj2PQQmdDr6SxVi2G6rQSF2qydo0NDmEiKTJ+eefX2e7T58+zJ49m6eeekrJmuTCCkJTnxeTHOsHfJbdcFo/j/6T1KQ6dMewDMchItu5NWvWsMsuu+Q6DNk+TQJuNrMyd3+rZqeZ9QNuBP6as8haKcNQq6rUpVqzJiJ5YORImDMnu9fs3x/GjGnaayZPnkxxcXGdfQ21R5s9ezbjxo3j97//fdMDFNl2twDHALPM7N/Ah8DuwCBgMWFgeEkj1ao1TUrJmpmd01gZd39i28MRkdZg8ODBjB07ts6+uXPn8u1vf3ursu+++y4nnHACI0eO5NRTT81WiCK13H2lmQ0EriYkbf0Jc2L/hDDw+xc5DK9Vqq6qZu2mMIOBOhk0LtWatXH17I9PjZWsiWRYU2u4cqVDhw7svffedfatXr16q3ILFizgyCOP5Ac/+AH33HNPlqIT2Zq7rybUsN2S41BavfLycio2VVDxZYVmMEhRqkN39E6yDABuAxYCX8tIdCLSar3zzjuUlpbyve99j5/97Ge5DkdEsqSsrAzaAJ1gU5dNdWY0kORS7WCwNMnupcDsaNy1q4Ez0hmYiLRe8+bN46ijjuLII4/khhtuYPny5bXHunXrlsPIZHtlZgcAFwD7AO0SDru7a6y1NNm1/66wETCoPrs6bEuD0tHB4BVCsiYJ1IBSJLk//elPrFixgokTJzJx4sQ6x9z1eyPZZWZfA14ClgB9gbeAXYCehM4G7+UsuFbos+LPwACDWJtY2JYGpfoYtCFfB9am4Tytlrony/Zk3LhxTJo0aav9AwYMwN3p1asXo0aNwt2TLiI5cBfwF+CrhDTiAnfvBRwNFAB35i601qe0V2lYcWhb2HbLttQr1d6gyRpcFgEHACcAD6UzKBERkSw6EDiXLZ3mCgDcfbqZ3Qncjdpmp01JjxKKviyiqKiIFy54gZIe6lzQmFQfg45Ksm8Tod3aTwg3soiISEvUBljn7tVmtgroHnfsXULFhKRRrDJGsRUrUUtRqh0M0vG4VEREJB8tIgyCC6G92vlmVvMs/zxgedJXiWRJSkmYmQ02s+J6jhWb2eD0hiUiIpI1k4DSaP0u4DhgDfA5YaSD0bkJSyRI9THoi0AJMCPJsX2i4wXpCkpERCRb3P3WuPWpZvZ14FSgAzDZ3V/IWXCtlGYwaJpUk7WGujO2BarSEIuIiEhWmVkb4HjgLXdfDODubwBv5DSwVqy8vJyKigoq1msGg1TVm6yZWS+gT9yuAUkehbYHzgf+m/7QREREMsvdN5vZH4FjCZO2S4bFz1hQUVFBWVmZkrVGNFSzdi5wK6ErswO/oG4Nm0fblcCITAUoIiKSYe8DXXMdxPaitLQ0fOJAUVFR2JYGNZSsjQPKCAnZdEJC9k5CmU3Af9x9VSaCExERyYL7gBvNbLq7f5rrYFq7kpISCp8sJLZTjDFPjVGtWgrqTdai+UCXApjZkcAsd9dMBSIi0tocBXQCFpvZ68DHUGe+QHf3c3MSWStU/kE5lbtUAjBy9kj6HdBP4601IqWhO9z9pXQkamZ2vZn928zWmNmnZvZMNHlufBkzs1Fm9pGZbTCzMjP7akKZtmb2CzNbaWbrzOxpM9sjocwuZjbezL6IlvFmtvO2vgcRSc2yZcu4+OKL2WOPPSgqKmL33Xfnoosu4sMPP6wtc/fddzNw4EB23HFHunTpwkknncTcuXNzGLVspw4HNgOfAntF299MWCRNypaUhRWDiqqKLdtSr5QHuzWzoWb2VzN7x8zeT1gWpXiaUuAR4DDCN5lKYKqZdYor8yPgGuAKYCCwAphiZjvElRlD6FZ9OuGXaEdgkpnFDx8yATiEMF7OsdH6+FTfr4g03+LFixkwYABz587l8ccf57333uN3v/sd8+bNY+DAgSxZsgQIDY2HDx/Oa6+9xvTp0yksLOToo49m1Sq1rJDscffejSx9Gj+LpCp+btAYMc0NmoJU5wY9HngGmArsC0wmjD/zDcKj0ldSOY+7D00479nAF9F5nrEw4/lI4B53fyoqcy4hYTsD+JWZ7QRcAJzn7lPizrOUMOnu82a2HyFBO9zdX4vKXAK8Ymb7uPu7qcQrIs0zYsQIYrEYU6dOpUOHDgD07NmTqVOn0rdvX0aMGMGzzz7L888/X+d148ePZ6edduLVV1/lpJNOykXosp0ws5eBi919Qdy+o4B/ufu63EXW+r099+3a9c2bN/P23Lf1GLQRqdas3Qw8TBiLBuAmdy8FvkoYDPe5Zl5/hyiGz6Pt3kA3oHYAQnffALxMqI0DOJQwj1t8mQ+A+XFlSoC1wGtx13oVWBdXRkQyYNWqVUyePJkRI0bUJmo1OnTowPDhw3nuuef4/PPPt3rtl19+SXV1Nbvssku2wpXt1+GEpzIARE9mphAGepcMemrWU2HFwlK7LfVKdVDcfYFbgGpCo8tCAHf/j5mNIiRzf2zG9R8E5gDl0Xa36OcnCeU+Ycu8bd0Ig/CuTFKmW1yZT929toGou7uZrYgrU4eZXQxcDKErsUg+Gjl5JHOWz8nqNft368+YY8ekXH7hwoW4O/vtt1/S4/vvvz/uzsKFCxk0aFCdY1deeSX9+/dX7zDJlYYGgJc0OfXQU3nh9ai+xcO2NCzVmrVqoDJKfj4FesYd+4jQILNJzGw04ZvNqe6eOAOCJxZPsm+rUyaUSVa+3vO4+1h3H+DuAwoLU81hRaQ+oVXD1mq+QyUev/rqq/nnP//JU089RUGBZq8Taa36HdCvNi1uU9QmbEuDUs1K3gV6ReszgZFm9iqhg8A1wJKmXNTMfgb8ADjS3d+PO7Q8+tkN+CBuf1e21LYtJzx67UxIHOPLvBxXpquZWU3tWtQergtb19qJtBhNqeHKlb59+2JmzJs3j1NOOWWr4/Pnz8fM2GuvLd/xrrrqKp588klefPFF+vRRW27JmmRf3hurGJBtVLakjK9/AOe8BbCZhZ2foOQy1aY3JNVk7fdAzTONWwkdDWr631cRGv+nxMweJCRqpfENOyOLCYnWMcC/o/LtCD0+r4vKzCJ0sT6G0OOTaNiO/djSRq0cKCa0XavZVwJ0pG47NhFJs06dOjF06FAeeeQRrrrqqjrt1tavX8/DDz/McccdR6dOoRP4lVdeyZNPPklZWRn77rtvrsKW7dNtZlbTpKamqvcOM0vsjqxx1tJox5fmUTYOiqrDdvUbY2FFN+jbN6dx5bOUkjV3fzhufZaZ9SP0tuwATHX3xJkNkjKzh4GzgVOAz82spv3YWndfG7UrG0MYSXoB8B/gJkJngQnR9b8ws18D90dt0D4DRgNvEZJI3H2+mU0m9B69iPBL+CtgknqCimTeQw89xGGHHcbRRx/NnXfeSd++fVm0aBE33ngj7s5DDz0EhF6j48eP529/+xu77LILy5eHyvXi4mKKixOnIpbtjZkNJ3xR7w7MA0a6e6OjD5hZX2A2YO5e3430X7ZUQtRYSug4l0i1bWm059MvUVS9JTsuqKqGUaNyGVLeazRZM7Mi4DJgmrvPBXD3D4HHmnG94dHPaQn7bwNGRev3ESaIfxjYBfgX8C13/zKu/FWER7ATo7LTgHMS2r6dCfycLb1GnwYub0bMItJEe+21FzNnzuT222/n7LPPZsWKFXTp0oXjjz+eiRMnssceYQzrRx55BIAhQ4bUef2tt97KKP3jvV0zs9MIndCGA/+Mfj5nZvu7+38beF0R8CShWcwR9ZVz915pDVhStqe3q113wA1swh/g4INzF1SupPg0odFkzd0rzOweYGhjZVM4V6M9baI2ZqPYkrwlK7ORMGjuFQ2UWQWc1eQgMyCuU6rIdqNHjx48+uijDZbR74Y04GpgnLvX3ERXmNmxhMqD6xt43b2EJy0v0UCytj0rL4eyMigthUx2vE68Ts32QZ9WEj910cs94ch3Tw+t4yWpVNuszQf6sKUBv4iISEZEtWOHAj9NOPQCDYyVaWYnACcSZqzJ6/EgNm3ag+9+N/3nracTdq3PPoOXX4aqKigogMGDYdddt+2c9V3npZe2XKd/f5gzJ2w/0rN9nbLzO8OQ3kP4Zs/tb1avUfXXS9WRarJ2C/Cgmc1y97cbLS0AmIbsERFJptDMZsZtj3X3sXHbnQm9/pONuXl0shOaWXfgUeA77v5lfUPH5Ivq6iIWJHax20apVFSvXBkSJgg/582Dzp237ZypXGfBgi3bT2wYySVcAkBFDJ7o14bRR96xXc5ikO5k7ceE3pVvmNkS4GMSxjRzd1U3i4hIKirdfUAK5Zoy5ubvgF+6++vbFFmWtG//PnPnZv+65eUwZAhUVEBREfztb5l5FJp4ndGjYeRI2LAB+PSA2v+J7jEKn3+Qkt9uf4laU6SarFUBKfX4FBER2UYrCX93EmeciR9zM9FRwBFmdmu0bUDMzCqB4Qk1d9utkhKYNi3zbdaSXadfv7C94ZYXiVWG/0GFDidv0vCnjUl16I7SDMchIiIC1HZsm0UYT/NPcYeOAeqbSDJxGPyTgRuBQcCytAfZgpWUZLZjQX3Xqdm+61cbYGmoIi2gmoodKjMfTAuneZVERCQfjQbGm9kM4FXgUmA34P8AzOxuYJC7DwGoGVqqhpkNAKoT9ydjZp2BDvFDgpjZJcABwPPuPik9b0kADt5zKSwNNWuVxMK2NCjVuUExs93NbLSZzTSzxWZ2QLR/pJl9LXMhiojI9sbdJwIjCQOjzyHMJX28u9f8Ze9OM+alrsdvgP+t2TCzm4FfEmbn+Xs05pukydudeuMYDlTShrc79c51SHkvpWTNzL4KvE2YfeAjwkTuRdHhPYErMxKdiIhst9z9EXfv5e5t3f1Qd3857tiwhga2dfdxDcxekGgAdQdrvxS4y913JQzQfnUzwpd6rFq1J+BRJwOPtqUhqdasPUAYa6038B2oMybFa8DX0xyXiIhItnQi6rgQPTXqBjweHfsbsE9uwmqdvrryz6H3B1BAFV9d+edch5T3Uk3WDgfucfe1bN1t+hO27rEjIiLSUnwG7BGtHwV85O4Lo+02NKHJkDRuSkW72segVRQydXOy6VglXqodDKobONYZ2JCGWESkFVm2bBm33XYb//jHP+rMDXrrrbfWzg06bNgwVq5cyaRJddtvz5w5k4EDB7J48WJ69eqVg+hlOzMVGBV1NLiGUJtWY1/CBO/SROXl5ZSVlVFaWkpJSQnl5eU88cQTLHzvbOBvtY9B/7PwO/z0p7DffrmNN5+lmqzNAM4Dnkly7PuEnjoiIgAsXryYww47jN69e/P444/Tt29fFi1axI033sjAgQMpLy9XEib55EeEQXXvBv4N3BZ37EzCRPKtRmISVR93p7Kyks2bN7N58+Y6640tb731FjfccAOVlZUUFBRw7LHH8txzz1FZWcn/sjtGNM4alZTyEtddV+8sYkLqydodwFQzewGYQHgUerSZXQl8GxicofhEpAUaMWIEsViMqVOn0qFDBwB69uzJ1KlT6du3LyNGjODZZ5/NcZQigbt/QhjDLZmjgY3pvuamTZs477zzqK6uxt2prq7OyvratWvpvHAhRxCqEJfuthsFBQVJE67Kym0b/+zrQClQVl3NM888U7t9Oo8BW8ZZ+4xdeewxOPDAbbpcizRoUGrlUh0U9yUzOwUYQ+jiDHAPsAQ4xd3/1dQARaR1WrVqFZMnT+bOO++sTdRqdOjQgeHDh3PzzTfz+eef5yhCkcaZ2f7AfkC5u69J9/l7VlVx0h//CBZmkTazuusJ+zDD4tfr2ZfstfHrhZ9/TgmhEV418K81a4h17Yq1bUvMjFgsRiwWw6KfMbM668mOJ9u3adkyur/zTu11Vu65J52XLiUGGEtreyk68HDxj2lzwcXp/ohblZQHxXX3Z4FnzWxvwpQfn7n7uxmLTES2NnIkzJmT3Wv27w9jxqRcfOHChbg7+9XTAGX//ffH3Vm4MLTfnjx5MsXFdUdYqK5uqJmsSHqZ2UNAobtfGm1/B5hImEx+jZkd4+7/Tuc1i834Trfs983bXFBAAdQ+hhwUi9HGLMzY7g7p+t37/HM87jrdVqyo3Y7vpWhAm7Wr03PNVqzJMxi4+3vAexmIRURakZpv+Incvc7xwYMHM3Zs3Wkb586dy7e//e3MBiiyxXHUbad2GzAJuIUwdNWtwInpvOCi9u1h0aJ0njIlbcrLqTryyNoZ1ttMnpyxmdwtmsndiopgzBgsmsl9q38Z9tQ4a41JOVkzs76EkaRLgN0Jc629BtwZJXAikmlNqOHKlb59+2JmzJs3j1NOOWWr4/Pnz8fM2GuvMPh8hw4d2HvvveuUWb16dRYiFanVjdCsBzPbA/gqcIG7v21mPwd+ncPY0qukhIIXX8ztTO5lZfDCC6HcnnvCkiWZiaEVSSlZM7NS4B+EITqeJYyt9hXgJOA0MzvW3V/KUIwi0oJ06tSJoUOH8sgjj3DVVVfVabe2fv16Hn74YY477jg6deqUwyhF6tgA1DyLPwJYA8yMttcCO+QiqIzJ9Uzu11+f+Wu3Mk2ZweANYE93P8fdr3P3c4BehDnbHshMeCLSEj300ENUVlZy9NFHM336dD744APKyso45phjcHceeuihXIcoEm82MCKavWAEMMXdaxpv9QY+zllkIqT+GHR/4LRoBoNa7v6lmd0L/CHtkYlIi7XXXnsxc+ZMbr/9ds4+++w6g+JOnDixdlBckTxxIzAZeBNYTZgbtMYphLFG02r9+vVuZts6oHwhsG3ja+S/1v4e26dSKNVk7UO2TNyeqIjQfk1EpFaPHj149NFHGywzbty4pPsHDBhQ2xFBJNPc/d9m1pMwW8HChKE6xgILk79ym665zVNYmdlMdx+Qjnjy1fbwHlORarJ2L3CbmZW7e21iZma7E3rJ3JWJ4ERERLLB3dcBs5Ls1+jNknOpJmtHEBpYLjKz19nSweDr0Xpp1AkBwN393DTHKSIiklFmdhCwD9Au8Zi7P5H9iESCVJO1w4EqQiPLPaMFtjS6/GZcWT27iOgxjohI/jOznQkjHZRA7ditUPfvWT4ma2MbL9LibQ/vsVGpTjfVO9OBtGb1DQ4qIiJ54S5gV0LFwyuEOa+/AM4nJHA/yF1o9XP3Vp/IbA/vMRXb3MBRRESkhRtKSNhej7Y/dPeyaIiqqcCVOYtMhCZON2VmPYAeJH+ePz1dQYmIiGRRd+B9d68ys43UHQT3L8CTuQlLJEh1BoM+wO+BQTW7op/x87IWpD06ERGRzFsO7BytLyU8+iyLtvdOUl4kq1J9DPoY0BMYCRwLHBktR8X9FBERaYn+SUjQAMYDt5rZr8zsYeB+4PlMXdjMOpnZX81snZktNbMz6il3gJk9b2YrzWyr3mtmdrmZzTSzTWY2LsnxIWa2wMzWm9mLZpaV2dOz8f7MrJeZuZmtjVtuztBbyolUH4MOBIa5+1OZDEZERCQHbgN2i9bvJ3Q2OA3oADwNXJHBaz8MVBCGw+oPPGtmb7r7vIRym4E/Ao8Af0tyno+AOwnt7+qMim9mnQmPcy8EngHuACYSht/KtIy/vzg7u3urnO2gKTMYVGQyEBERkVxw90XAomh9M3BNtGSUmXUETgUOiKZz/KeZPQ2cDfxvQozvAu+aWdLHsu7+l+icA4DE+dy+A8xz9z9FZUYBK81sX3dfkMa3VEcW31+rl+pj0LuAH0cfvIiIiGy7/wGq3P0/cfveBL6a5ut8NTovUDtbw6IMXCdRtt5fjaVm9qGZ/TaqTWw1Uh1nbbyZ7QssiWYw+HzrIpq1oD6ucYJFRPKKmZ3flPLu/psMhFFMGM8t3hfU7Y2arut8moXrJLtuNt7fSkJzrTmER9gPEzpFDk3zdXIm1d6gw4DrCbMYHMLWj0SVjYiISEvyGFv+djU2crkDmUjW1gI7JuzbEfiyhV4nJ9eNHrHOjDY/MbPLgY/NbEd3X5POa+VKqm3WbgP+Clzg7qszF07rYo3+/ouISA6tBf5M6AG6OAfX/w9QaGZ93X1htO8gILHx/baaB9Q+/YqaNO2Vgeskytb7S5RqEt5ipNpmbVfgESVqIiLSSvQGfkqYYmoqYe7PIcAqd1+auGQigKjt2F+A282so5l9AziZkDzWYUE7oCjabmdmbeOOF0bHC4CC6HhNhcxfgQPM7NSozC3AW5nsXJDN92dmXzOzfcwsZma7Aj8Hytw98RFsi5VqsvZPYL90XNDMBpvZ02a2LBoXZVjC8XHR/vjl9YQybc3sF9F4LOui8+2RUGYXMxtvZl9Ey3gLk/WKiMh2LkrC7nD3/wEGA/MJw3YsN7M/mNlxZpaNKRmHE4aiWAH8AbjM3eeZWc9ovLCeUbk9gQ1sqZXaALwbd56bon3/C5wVrd8E4O6fEnpl/oTQ5vxrZG++04y/P6APMJnweHUusAk4PWPvKAdSfQx6JfBHM/uc8IEkdjDA3atTPFcx4cN8IlqSmUro2lsjsY3cGEJ2fjrwGTAamGRmh7p7VVRmAmEg3+MIVaKPEbL5k1KMU0REtgPu/hrwmpn9kPA34lzC+Gp/ApIO4prGa68CTkmy/7+Ev5c120to4LGeu48CRjVwfCqwb7MDbaZsvD93/wMhEWy1Uk3W5kc/60uuPNVzufs/gH9AqEWrp9gmd1+e7ICZ7QRcAJzn7lOifWcTpgg5GnjezPYjzLRwePRLiJldArxiZvtE47mIiIjE2xXoRajlKSD0MhTJuVSTtdvJbo/Pw81sBbAaeAm40d1XRMcOBdoAL9QUdvcPzGw+cBhhWpASQsPR1+LO+SqwLiqzVbJmZhcDFwMUFRWl+e2IiEg+MrP2hEFjzyZ84f+QMOzD9/XFXvJFqrVhozIcR7zJhAaJiwnfcO4EpkePODcB3QhDiCR+4/kkOkb081N3r00w3d2jBLAbSbj7WGAsQMeOHTUUiYhIK2ZmRxMStG8TKiP+Ahzj7i/mNDCRJFKtWatlZsWEquKPomk50srdn4zbfNvMZhEecZ5A+GWqNzTq1v4lS7gSy4iIyPbpBWANYeiOvwDrCZ0Sj0pW2N2nZzE2kTpSTtbM7ETC49CDol0Dgdlm9hgw3d0nZCA+3P0jM/sQ6BvtWk5oS9CZuiMydwVejivT1cyspnbNzAzoQqiBExER2REYRtwYZNRt5O5s+ZJfkL2wROpKdQaDU4CngGnAj4H74g4vJtzoGUnWovm9dgc+jnbNAjYDx9RcMxq2Yz+2tFErJ/QyKYnbVwJ0pG47NhER2T4dmesARFKVas3arcBv3f3CaBC6+GRtLmEclZREj1H3jjZjQE8z6w+sipZRhMTwY0KbtbsJ47P8FcDdvzCzXwP3R23QaobueIsw5AfuPt/MJgO/MrOLCN+MfgVMUoNREZGWwcyGA9cB3Qnjb41091fqKVsKXAUMAnYC3gPG1Denp7u/lIGQRTIi1QH/9gMmRuuJbb4+J7RhS9UA4I1oaU+YyuoNwiPWKqAf8HfCNBWPE3pulrh7/FxiVxHaGEwk9PJcC5wUN8YawJnAm4R2Cc9H6/Fjt4mISJ4ys9OAB4G7gIMJT0WeixtENdFhwNvAd4EDgF8CY80so+OkiWRDqjVrawhtxJLpRd22Yw1y9zIanq9raArn2AhcES31lVlFGOVYRERanquBce7+aLR9hZkdC1wGXJ9Y2N3vStj1SzM7kjByf0aa6YhkS6o1a1OA6xOma/Jo3q7LgefSHZiIiGyfzKyIMKbmCwmHXiDUoKVqR5LMuCPS0tSbrJnZ+2ZW0/PzRsL4ZO8Spm1ywvxcc4A9aGCKCxERkQSFZjYzbrk44XhnQu/LxN778eNpNigawWAI0fiZIi1ZQ49BewFtIczZZWaHENqXDSW0LRtMGMD2Fnf/KMNxiohI61Hp7gNSKJfYRjqlsTLN7BuER58/dPcZzYhPJK+kPM6au39ImJNTREQkk1YSKgUSa9G60shYmWZ2OGH+6Vvc/ZeZCU8kuxprs6bR/kVEJKvcvYIwpuYxCYeOoYGxMs1sMKEN9W3uPiZjAYpkWWM1a7eZWeIcnMm4u5/beDEREZGUjAbGm9kMwhBNlwK7Af8HYGZ3A4PcfUi0XQo8CzwC/N7Mamrlqtw95RELRPJRY8laf2BTCudRDZyIiKSNu080s12BmwiD4s4Fjnf3pVGR7sBecS8ZBnQAro2WGksJbbBFWqzGkrVT1DhTRERywd0fIdSUJTs2LMn2sGRlRVq6VMdZExEREZEcULImIiIikseUrImIiIjksXrbrLm7EjkRERGRHFNCJiIiIpLHlKyJiIiI5LGUp5uSpuuy7FM+uQ92ued8KEicp1hERESkcUrWMmjXFavpuh4WHz2A3gd9PdfhiIjkhwceyHUEIi2KkrUsWPDd4+h9yY25DkNEJD8oWRNpErVZExEREcljqlnLgkVfvsf0xdNzHYaIiIi0QErWMqgoVgTAbxeOY/YT43IbjIiIiLRIStYyqGe7ngDcfNBNdDrhmBxHIyKSH44YdUSuQxBpUZSsJbFhQ28OPhjMtiyxWN3txvbFYnDIsr6UAP+zw37sv+fgXL8tERERaYGUrCURi22mRw9wr7tUVze+r7Iy7KuqgvnzdwVg2Udt2T/H70lERERaJiVrSbRt+yFPP73t53n81HfhL+DV234uERER2T5p6I5sMMt1BCIiItJCKVkTERERyWNK1kRERETymJK1rNBjUBEREWkeJWsZpKZqIiIisq2UrImIiIjkMSVrIiIiInlMyVomuec6AhEREWnhlKxlVEjWql2N10RERKR5lKxlgamngYiIiDSTkjURERGRPJb1ZM3MBpvZ02a2zMzczIYlHDczG2VmH5nZBjMrM7OvJpRpa2a/MLOVZrYuOt8eCWV2MbPxZvZFtIw3s50z/w5FRERE0icXNWvFwFzgSmBDkuM/Aq4BrgAGAiuAKWa2Q1yZMcCpwOnAN4EdgUlmVhBXZgJwCHAccGy0Pj6db0REREQk0wqzfUF3/wfwDwAzGxd/zELjrpHAPe7+VLTvXELCdgbwKzPbCbgAOM/dp0RlzgaWAkcDz5vZfoQE7XB3fy0qcwnwipnt4+7vZvp9xlOnUBEREWmufGuz1hvoBrxQs8PdNwAvA4dFuw4F2iSU+QCYH1emBFgLvBZ37leBdXFlRERERPJeviVr3aKfnyTs/yTuWDegCljZSJlP3bfUaUXrK+LK1GFmF5vZTDObWVlZ2fx3EH9OVKUmIiIi2ybfkrUaiVmOJdmXKLFMsvL1nsfdx7r7AHcfUFiY9afDIiIiIknlW7K2PPqZWPvVlS21bcuBAqBzI2W6WtwAZ9F6F7autcs8DbMmIiIizZRvVUiLCYnWMcC/AcysHaHH53VRmVnA5qjMhKjMHsB+bGmjVk7odVoSt68E6EjddmxZMXr0blz9RLavKiIiIq1B1pM1MysG9o42Y0BPM+sPrHL3/5rZGOBGM1sA/Ae4idBZYAKAu39hZr8G7jezFcBnwGjgLWBqVGa+mU0m9B69iFC39StgUjZ7gnbv/gUAe+21kZ17ZOuqIiL5bf78XEcg0rLkomZtAPBi3PZt0fI4MAy4D2gPPAzsAvwL+Ja7fxn3mquASmBiVHYacI67V8WVORP4OVt6jT4NXJ7m99KgzZtXAXDrrevodmw2rywikr80A59I05hrELCtdOzY0detW7dN51i9ejV3DBjAA4sWUTVzJgWHHpqm6EREWjYzW+/uHVMoN5zQBKY7MA8Y6e6vNFC+H/AQMAhYRXiicofrD520cPnWwaBFq66uZvbs2dxxxx0MHDiQJYsXA1BQUNDIK0VEJJ6ZnQY8CNwFHExob/ycmfWsp/yOwBRCJ7KBwA8Jid7VWQlYJINUs5bEboWFPuHEE3F3qKoCd7y6Gqqra39WVVZSUVHB5o0bWbd2LWvXrGHNF1+wadMmCoDdunfn+L33pusrr8Abb0D//rl+WyIieSGVmjUz+xfwlrtfFLdvIfBnd78+SfnLgHuBr0SDqWNmNwGXAXuodk1aMiVrSQww85nNfG21GRaLYWYQi0GnTvDmm9C1a1pjFBFpqcysAng7btdYdx8bd7wIWA+c7u5/itv/MHCAux+R5JxPALu6+wlx+wYCM4A+7r44/e9EJDvybeiOvPBuu3Z8MX8+VlBArGYpLKzz02IxrKAgtJSNxWpbzOq5sohIoyrdfUADxzsTxtNMNpvN0fW8phvwYZLyNceUrEmLpWQtiepYjJ169cp1GCIi27umzmaTrHyy/SItiiqCREQk36wkzAHd0Gw2iZbXU54GXiPSIihZExGRvOLuFYTZao5JOHQM9c9CUw58M5r1Jr78R8CSdMcokk1K1kREJB+NBoaZ2YVmtp+ZPQjsBvwfgJndbWbT4spPIHRKGGdmB5jZd4D/BUarJ6i0dGqzJiIiecfdJ5rZroQpB7sDc4Hj3X1pVKQ7sFdc+S/M7BjC7Dczgc+BBwhJn0iLpqE7kkjHDAYiIpJcqjMYiEigx6AiIiIieUzJmoiIiEgeU7ImIiIiksfUZi0JM6sGNuQ6jhQVApW5DiJFijVzWlK8ijUzWlKs7d1dlQUiKVJv0ORmNzIVSt4ws5mKNf1aUqzQsuJVrJnR0mLNdQwiLYm+2YiIiIjkMSVrIiIiInlMyVpyY3MdQBMo1sxoSbFCy4pXsWaGYhVppdTBQERERCSPqWZNREREJI8pWRMRERHJY0rWRERERPJYq0/WzGy4mS02s41mNsvMvtlI+X5m9pKZbTCzZWZ2i5lZQpkjonNtNLP3zezSbMdqZqVm9ncz+9jM1pvZW2Z2fpIynmTZN8ux9qonjmMTymXkc21GvKPqidfNrGtUJiOfrZkNNrOno3vPzWxYCq/JyT3b1Fhzec82I9ac3bPNiDWX9+v1ZvZvM1tjZp+a2TNmdkAKr8vZv7MiLVGrTtbM7DTgQeAu4GDgNeA5M+tZT/kdgSnAJ8BA4IfAdcDVcWV6A/+IznUwcDfwCzM7NZuxAocBbwPfBQ4AfgmMNbMzkpT9KtA9blmY5VhrHJsQx/S4c2bkc21mvD9NiLM78BJQ5u4rEsqm9bMFioG5wJWkMItGLu/ZpsZKDu/ZZsRaIxf3bFNjzeX9Wgo8Qvh/exRhBoWpZtapvhfk+J4VaZncvdUuwL+ARxP2LQTurqf8ZcAawlQoNftuApaxpefsvcDChNc9BpRnM9Z6zvFH4Km47VLAgc45/lx7RXEMaOCcGflc0/HZAj2AKuCMTH+2CdddCwxrpEzO7tmmxlrP67Jyzzbjc83pPbstn2uu7tfoOsXRtU9qoExe3LNatLSkpdXWrJlZEXAo8ELCoRcI3wKTKQFecff4b7PPA7sR/vGuKZN4zueBAWbWJouxJrMj8HmS/TOjR0/TzOzI5sRYYxtj/YuZrTCzV83suwnH0v65piHeGhcAq4GnkhxL22fbTDm5Z9Mo4/fsNsr6PZsGubxfdyA8sUn2/7RGS79nRbKu1SZrQGeggFDVHu8ToFs9r+lWT/maYw2VKYyu2RzNibUOMzsRGELdwSY/JnyLPRX4DvAuMM3MBjczzubGuha4Fvg+cDwwDZhoZmfFlcnE59rceGuZWQw4H3jC3TfFHcrEZ9scubpnt1kW79nmyOU922x5cL8+CMwByhso02LvWZFc2R4mck8c9deS7GusfOL+VMo0R1NjDYXMvgFMAH7o7jNqT+b+LuEf5RrlZtaL8Efo5WzF6u4rgQfids00s87Aj4DfNXLOZPubo1mfLXAc4bHSY3VOltnPtqlyec82S47u2ZTlyT3bHDm7X81sNHA4cLi7VzVSvMXdsyK51Jpr1lYS2k4k1p50ZetvbDWW11OeuNfUV6YS+KxZkTYvVgDM7HDgOeAWd/9lCtf6F9C3OUFGmh1rI3Fk4nOFbY/3YuA1d5+XQtlt/WybI1f3bLPl4J5Nl2zds9siJ/ermf0MOB04yt3fb6R4i7tnRXKt1SZr7l4BzAKOSTh0DKGHUTLlwDfNrF1C+Y+AJXFljk5yzpnuvjmLsRI9wngOuM3dx6R4uf6ERyLN0txYU4gj7Z8rbFu8ZrYbcALwaIqX6882fLbNlJN7trlycc+mUX+ycM82V67uVzN7EDiDkKgtSOElLeqeFckLue7hkMkFOA2oAC4E9iO0p1gL7BkdvxuYFld+J8I3uicJQwt8h9Br6Zq4Mr2BdcCY6JwXRtc4NcuxlkZx3E/4BlqzdIkrMxI4hfDt+avRORz4TpZjPZfwj/l+wD6ERy8VwFWZ/lybE2/c624CvgA6JDmWqc+2mPBHtD+wHrglWu+Zh/dsU2PN5T3b1Fhzds82NdYc368PR/fbUQn/T4vjyuTNPatFS0tdch5Axt8gDCd8W9tEqGEZHHdsHLAkoXw/QhuOjYRvnbcSdSePK3MEMDs652Lg0mzHGm17kiW+zI+A9whjNa0CXgGOz0Gs5wLvRP/4rgFmAmclOWdGPtdm3gcWxfBIPefLyGfLliEWEpdx+XbPNjXWXN6zzYg1Z/dsM++BXN2vyeJ0YFQjv185+3dWi5aWuNSMaSMiIiIieajVtlkTERERaQ2UrImIiIjkMSVrIiIiInlMyZqIiIhIHlOyJiIiIpLHlKyJiIiI5DElayLbyMw8hWWJmfWK1oflQcy9EuIrbcJrb4p73YeZi1JERGD7mMhdJNNKErb/CrwJjIrbt4kw+GcJsCg7YaXkTuBZwgCwqfotMBW4GTgoE0GJiMgWStZEtpG7vx6/bWabgJWJ+yPJ9uXSonrirJe7LwOWmdmnGYpJRETi6DGoSJYkewxqZuPM7EMzG2Bmr5nZBjN718xOiI5fHT1CXWNmfzezLgnnLDSz681sgZltMrOPzOyBhEmymxrn0CiWL8xsbRTPLc1+4yIisk1UsyaSezsCTwA/BT4CbgSeMrOHgf8BRgBfIUxq/TDw/bjX/g44CbgXeI0w6fUdQC/g1KYGYmZ9gKeBPwO3EybP7gv0afK7EhGRtFCyJpJ7OxAmqX4ZwMw+IrR5OxHY392rov0HAFeYWYG7V5nZN4HTgHPd/YnoXFPNbBXwOzPr7+5zmhjLIUARcJm7r4n2Td+WNyciIttGj0FFcm9dTaIWWRD9nFqTqMXtLwS6R9vHEmq+nooehxaaWSHwQnR8cDNimQNsBp40s++aWddmnENERNJIyZpI7q2O33D3imj184RyNftr2qN1JdSCrSUkWDXLiuj4rk0NxN3fA4YS/m0YDyw3s3+Z2RFNPZeIiKSHHoOKtFyfARuBb9Zz/KPmnNTdXwReNLO2wDcIbdeeNbNe7r6yWZGKiEizKVkTabkmAz8GdnL3aek+ubtvAqabWTHwd6A3oGRNRCTLlKyJtFDuXmZmfwD+bGajgRlANaEn6PHAj939P005p5ldSmjr9g/gA6AzcD2hlm5u+qIXEZFUKVkTadnOAq4AzicM+bEJWAI8D3zSjPO9CRwH3E1oE7cK+CdwprtvSEO8IiLSRObuuY5BRLLMzHoBi4ELCGO8VXmK/xiYmQEFwK+BIe6+R6biFBER9QYV2d79mtCDtCm9PW+MXnNORiISEZE6VLMmsh0ysyLgwLhd77r7lym+tjuwe7RZ4e5vpTs+ERHZQsmaiIiISB7TY1ARERGRPKZkTURERCSPKVkTERERyWNK1kRERETymJI1ERERkTz2/wFhVd2tOFlaOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "info={}\n",
    "info['label1'] = {'label':'Temperature','units':' [K]'}\n",
    "info['label2'] = {'label': 'H2','units':' Mass Fraction'}\n",
    "info['label3'] = {'label': 'O2','units':'  Mass Fraction'}\n",
    "info['label4'] = {'label': 'OH','units':'  Mass Fraction'}\n",
    "info['loc_x'] = 0.45\n",
    "info['loc_y'] = 0.75\n",
    "info['xlim'] = [0,1]\n",
    "\n",
    "info['inset_x1'] = 0.55 + 0.5 \n",
    "info['inset_y1'] = 0.4\n",
    "info['inset_x2'] = 0.3 \n",
    "info['inset_y2'] = 0.4\n",
    "\n",
    "spNumber = 53\n",
    "info['xlim2'] = [IgnTime[spNumber]*0.95,IgnTime[spNumber]*1.05]\n",
    "info['xlim'] = [0,time_simulation[-1,spNumber]]\n",
    "\n",
    "x  = time_simulation[:,spNumber]\n",
    "y1 = solution[:,spNumber, T_indx]\n",
    "y2 = solution[:,spNumber, H2_indx]\n",
    "y3 = solution[:,spNumber, O2_indx]\n",
    "y4 = solution[:,spNumber, OH_indx]\n",
    "helper.makefigure(x, y1, y2, y3, y4, info, fs_label=16, fs_tick=14)\n",
    "plt.savefig('H2_TimeProfile.pdf', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finalize Kokkos. This deletes the TChem object and the data stored as Kokkos views"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "gkde = stats.gaussian_kde(dataset=IgnTime)\n",
    "x = np.linspace(start=np.min(IgnTime),\n",
    "                stop=np.max(IgnTime), num=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAEGCAYAAAAZjzycAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAl4UlEQVR4nO3dd3hUZd7/8fc3lVQQEnoHpYQiEBBEBBQbuopYFpUV9rfo4qKiYleUx4KsdWERlEeUtTzqguiKshYEDE0pUgKCLmGlCAoGSUibkMn9+2OGEEbKhMzJPeX7uq5cZMo59ydD+HBmzjn3EWMMSimljoiyHUAppYKNFqNSSvnQYlRKKR9ajEop5UOLUSmlfMTYDlBZWlqaadmype0YSqkws2bNml+MMen+Pj+oirFly5asXr3adgylVJgRke1Veb6+lVZKKR9ajEop5UOLUSmlfGgxKqWUDy1GpZTyocWolFI+tBiVUsqHFqNSSvnQYlRKKR9BdeaLCh5Dhix1dP0ffHCOo+tXqjp0i1EppXxoMSqllA8tRqWU8qHFqJRSPrQYlVLKhxajUkr50GJUSikfWoxKKeVDi1EppXxoMSqllA8tRqWU8qHFqJRSPrQYlVLKhxajUkr50GJUSikfWoxKKeXD0WIUkTtFZJOIbBSRt0WklpPjKaVUIDhWjCLSBLgdyDTGdAKigWFOjaeUUoHi9FvpGCBBRGKARGC3w+MppVS1OVaMxpgfgWeBHcAeIM8Y85nv80TkZhFZLSKr9+3b51QcpZTym5NvpU8DrgBaAY2BJBEZ7vs8Y8wMY0ymMSYzPT3dqThKKeU3J99KDwL+a4zZZ4w5BMwFznZwPKWUCggni3EH0FtEEkVEgPOBzQ6Op5RSAeHkZ4xfA3OAb4Bs71gznBpPKaUCJcbJlRtjHgUedXIMpZQKND3zRSmlfGgxKqWUDy1GpZTyocWolFI+tBiVUsqHFqNSSvnQYlRKKR9ajEop5UOLUSmlfGgxKqWUDy1GpZTyocWolFI+tBiVUsqHFqNSSvnQYlRKKR9ajEop5UOLUSmlfGgxKqWUDy1GpZTy4eg1X1R4MMZQUrKXwsI9xMQkkJDQgPj4OrZjKeUYLUZ1TOXlh9izZwnbt3/Er79+S1lZYcVjbdveQEbGaADcbhcgREfHWUqqVOBpMarfKC4uZsGC6yku/qnivri42iQlNcXtdpGU1Kji/t27v2TLllfIyLiFRo0G4LmEuFKhTYtR/UZCQgItW17Orl0LaNnyCho3HkCtWnWP+dz9+7MpKtrDqlWPUK/emfTo8QgJCek1nFipwNKdLwpjDFOnTmXevHkV951++vUMHDiL1q2HHrcUATp3HkuXLuOIi6tNbu46vvxyFLm5G2oitlKO0WKMcKWlpYwePZrbbruNG264gb179wIgEu3X2+KoqBhatRrCeee9SVpad1yu/SxbNpYdOz5xOrpSjtFijGAlJSVcfvnlzJgxg/j4eF566SXq169/SuuKj69Dnz7P0br1NRhTxqZNL3LoUEGAEytVM/QzxghVUlLClVdeyaeffkr9+vWZN28evXr1qtY6o6Ji6Nz5dpKSGpOW1p3Y2OQApVWqZmkxRiCXy8XQoUP55JNPSE9PZ+HChWRkZARs/a1bX33U7UOHComNTQrY+pVymr6VjkBr1qxh4cKFpKWl8cUXXwS0FCszxpCTM4cFC4ZRULDLkTGUcoIWYwQ6++yzWbRoEZ999hmdO3d2cCTDL7+sobT0ACtXPsChQ0UOjqVU4GgxRpD8/PyK7/v06UO3bt0cHU8kiu7dHyY5uQUHD/7A2rVPYoxxdEylAkGLMUKsX7+eVq1aMXfu3BodNzY2ibPOeoqYmGT27Mli+/YPa3R8pU6FFmMEyMvL4+qrr2b//v3Mnz+/xsdPTm5G167jANi48e8cPPhDjWdQqiq0GMOcMYY//elPbN26la5du/L3v//dSo6mTQfRtOlFuN0u1q9/zkoGpfylh+uEuRdffJH33nuPlJQUZs+eTUJCgrUsXbrciTFuOnS4yVoGpfyhxRjGvv32W+655x4AZs6cyemnn241T2xsEpmZj1rNoJQ/9K10mDr8FrqkpISRI0dyzTXX2I50FGMMs2fPpqyszHYUpX5DizFMiQjTpk1j8ODBTJ482Xac37jxxhu59tprgzKbUo4Wo4jUEZE5IrJFRDaLSB8nx1NH69atGx9//DGpqam2o/zGsGHDABg/fjzbt2+3nEapozm9xTgZ+MQY0x7oCmx2eLyI53K5+PTTT23HOKlLL72U3//+9xQXF3PXXXfZjqPUURwrRhFJBc4FZgIYY0qNMQecGk95PPnkk1x88cU88MADtqOc1HPPPUdSUhJz587lk090/kYVPJzcYmwN7ANeE5G1IvKKiPxmihURuVlEVovI6n379jkYJ/ytX7+ep556CoDBgwdbTnNyTZo0YcKECQDceuutlJSU2A2klJeTxRgDdAemG2O6AYXA/b5PMsbMMMZkGmMy09P1WiGnyu12M2rUKMrKyhgzZgz9+vWzHckvY8eOpWPHjuTk5PDyyy/bjqMU4OxxjLuAXcaYr72353CMYlT+GTJk6Qkf37btPbKzV1OrVn3++9/LT/p8Xx98cE514p2y2NhYpk2bxsqVKxk9erSVDEr5cqwYjTE/ichOEWlnjPkOOB/41qnxIllx8S9s3jwD8FycKjY20XKiqunfvz/9+/e3HUOpCk6f+XIb8JaIxAHbgD86PF5E+v77WZSVFdGwYV8aNQqNt9DHk5uby/79+62fpaMim6PFaIxZB2Q6OYaCjh3/QnR0Aq1bXxXSF7xfvnw5l112GWeccQbLly8nKkrPP1B26G9eGIiNTaRTpzEkJja0HaVaOnfuTHx8PF9//TXvvPOO7TgqgmkxhrDc3PWUlYXPIS4pKSk8+eSTANx///0UFxdbTqQilRZjiCou3suKFfewcOFwSkvzbMcJmBEjRtC1a1d27tzJ3/72N9txVITSYgxRmzZNw+0upk6d9sTF1bYdJ2Cio6N57jnPRLYTJ07k559/tpxIRSItxhCUm7ueH3/8gqioODp1utV2nIA7//zz+d3vfkdBQQHjx4+3HUdFIJ2oNsQY4yY72zNV1+mn3xDyO1yO5+mnn2bfvn0MHz7cdhQVgbQYQ8yOHf8mL+8/1KpVn7Ztr7cdxzHt27dnxYoVtmOoCKVvpUPIoUNFFWe4ZGTcQkxMLcuJak5hYaHtCCqCaDGGkJiYBLp0uZumTS+kSZPzbcepEbm5uQwbNow+ffroZRBUjdFiDCEiQuPG59Kjx/iQPsOlKpKSkvjqq6/Izs7mtddesx1HRQgtxhBRXPyL7QhW1KpVi7/+9a8APPzwwxw8eNByIhUJtBhDwNKlS/n886vZtGma7ShWXHvttfTu3Zu9e/cyadIk23FUBNBiDHLl5eXcddddGOMmOjredhwrRIQXXngB8FwOQS+epZymxRjk/vnPf7Jq1Sri4+vStu11tuNY07t3b4YNG4bL5eL++3W+Y+UsLcYgVlJSUlECHTrcRExMaE1AG2iTJk0iKSmJunXrUl5ebjuOCmN6gHcQmzp1Ktu3b6dTp040b36J7TjWtWjRgu3bt1OvXj3bUVSY0y3GIJWbm8sTTzwBwLPPPotItOVEwUFLUdUELcYg5XK5uOCCC7jwwgu56KKLbMcJKsYY5s2bx7nnnkt+fr7tOCoMaTEGqcaNGzN79mw+/PBD21GC0sSJE1myZEnFVrVSgXTCYhSRa7x/tqqZOAo8W0SHxcdH5iE6JyIiTJkypeIwns2bN9uOpMLMybYYH/D++Z7TQZTHqlWr6Nu3L0uXVu260JGmZ8+ejBo1irKyMm677baj/jNRqrpOVoy5IrIIaCUiH/p+1UTASGKM4d5772XFihXMmzfPdpygN3HiROrWrcsXX3zBnDlzbMdRYeRkxXgpnq3GX4DnjvGlAmj+/PksXryYunXr8sADD5x8gQiXlpZWcfGsO++8U8+jVgFzwmI0xpQaY74CzjbGfOn7VUMZI4Lb7a44mPuhhx6iTp06dgOFiJtuuonMTM+ly3NyciynUeHihAd4i8g8wHi//83jxpjLnYkVed544w02btxIixYtGDNmjO04ISM6Opp//vOfpKWlkZKSYjuOChMnO/PlWe+fQ4GGwJve29cBPziUKeIUFxdXXPTpiSee0D3RVdSqlR40oQLrhMV4+O2yiDxujDm30kPzRCTL0WQRZP369eTn53PmmWdy/fXhex0XpxUWFvLYY4/RrFkzbr01/K6eqGqOv+dKp4tIa2PMNqg4rjHduViRpXfv3mzbto2ffvqJqCg95v5UZWVl8fTTT5OQkMAll1xCmzZtbEdSIcrff4V3AotFZLH38J1FwB2OpYpA9erVIyMjw3aMkHbJJZdw/fXXU1xczKhRo3QGHnXK/C3GxcDLwK94dsa8DOhe6WrauXMnU6ZMweVy2Y4SNiZPnkx6ejqLFy9mxowZtuOoEOVvMb4OtAKmAI97v3/DqVCRYsKECYwdO5a7777bdpSwkZaWxosvvgjAPffcw7Zt2ywnUqHI32JsZ4wZZYxZ5P26GTjDyWDhbtOmTcyaNYuYmBjGjh1rO05Yueaaa7jmmmsoKCjgxhtvxO12246kQoy/xbhWRHofviEiZwHLnIkUGR588EHKy8u56aabaNu2re04YWf69Ok0atSI9PR0ioqKbMdRIcbfvdJnATeKyA7v7ebAZhHJBowxposj6cLU0qVL+fDDD0lKSuKRRx6xHScs1atXj1WrVtG4ceOIuQa3Chx/i/FiR1NEEGMM9913HwDjxo2jYcOGlhOFryZNmlR873K5KCsrIykpyWIiFSr8KkZjjF6vMkA+//xzli9fTnp6OuPGjbMdJyLk5ORw7bXX0rFjR15//XXdglQnpRfDqmGDBg3ijTc8O/RTU1Mtp4kMpaWlbNmyhW+++Ybzzz+fkSNH2o6kgpzjp1mISLSIrBWRj5weKxRERUUxfPhwhg8fbjtKxOjQoUPFITxjxowhOzvbciIV7Gri/LOxQMTPPV9UVMSOHTtO/kTliBEjRnDjjTdSVFTE0KFDOXDggO1IKog5Wowi0hTPZLevODlOKJgyZQpnnHEGU6dOtR0lIokI06dP58wzz2Tr1q0MHz5cTxlUx+X0FuPfgHuBiP4NzM3NZdKkSbhcLtq3b287TsRKTExk7ty51K1bl48//livwKiOy7GdLyJyGbDXGLNGRAac4Hk3AzcDNG/e3Kk4Vk2cOJG8vDwGDRrEoEGDbMeJaK1ateLtt99mw4YNXHHFFbbjqCDl5F7pvsDlIjIYqAWkisibxpij9joYY2YAMwAyMzND5lJvQ4b4dxW/oqI9fPHF3wE4dGiY38t98ME5p5xNndiFF17IhRdeaDuGCmKOvZU2xjxgjGlqjGkJDAMW+pZiJNi8+RXKyw/RtOkF1KnTznYc5WPbtm0MGDBArxejjqKzojooL+8/7Nr1OVFRsbRvP8p2HHUMDz/8MF9++SUXX3wxe/futR1HBYkaKUZjzGJjzGU1MVYwiY6Op0GDs2jZcghJSY1tx1HH8NJLL9G9e3e2bt3K4MGD9RKsCtAtRkclJzend+9nyMj4i+0o6jhSU1OZP38+bdq0Yc2aNVx55ZWUlJTYjqUs02J0gDEGY47sR4qK0jMvg1mDBg349NNPadCgAV988QVXX301paWltmMpi7QYHbBr16csX34HBw58bzuK8lObNm1YsGAB9erV02MclU4iEWhut4vNm/+X4uK95Odvo04dneg8GB3vsKlOnZ5m//4NvPlmQ958079Dq45FD7cKbVqMAbZt2xyKi/eSmtqGZs0usB1HVVGdOmcc9Z9ZUdFPxMfXJTo6zmIqVdP0rXQAuVwH+P57z5RiGRljEIm2nEhVR2HhbpYsGcPKlQ9QVlZsO46qQVqMAfTdd69RVlZI/fpnUb9+T9txVDWVlRVRXl7K3r0rWb78DkpL82xHUjVEizFADh7cwQ8//AuI0sNzwkTt2m3p1+9FEhIa8uuv37JkyRiKi3+2HUvVAC3GAMnP30pUVCwtWgwmNbW17TgqQJKTm9Ov3zRSUlpTULCdL7/8MwcOfGc7lnKYFmOANGlyHuef/zYdOtxsO4oKsISEdM45Zyppad1wuXJZuvQ2Skr2246lHKR7pQMoISHNdgTlkLi4FPr0eY61a/9KamoratWqazuScpAWYzX9+OMiDh06SIsWl+pe6Grwdzq2UxWI4wqjomLp3v2ho+7Lz88hMbEJMTG1qr1+FTy0GKuhrKyY7OzJuFy5xMXVoXHjc21HUg6rfOnVoqKfWLbsDmrVSqNnz8dJTm5qMZkKJP2MsRr+85+3cLlyqVOnA40a6ZkOkcbtdhEbm0x+/la+/HIUe/Zk2Y6kAkSL8RQVFf3E1q1vA9C5822I6EsZaVJSWtC//ys0atSfsrJCVq58iOzsKbjdOgFFqNN/zado06bplJeX0qTJIOrW7Ww7jrIkNjaJnj0frzjTadu22SxZMprvv9cJREKZFuMpyMrKYvfuhURHx5ORMdp2HGWZiNC27TD69ZtGYmJj8vJy2Ldvn+1Yqhp058spmDFjBgBt295AQkIDy2lUsDjttI4MGDCTfftW07dv34r7CwoKSE5OtphMVZVuMZ6CWbNm0bXrPbRte53tKCrIxMYm07jxgIrb//73v2ndujVz5861F0pVmW4xnoKYmBhatrzcdgwVAv7v//6Pffv2cdVVV9G06QV07jyWuLjaAR9H538MLN1irILXX3+dPXv22I6hQsg//vEPpk6dSnR0PLt2fc7ChTeye7ce1hPstBj9tHbtWkaOHElGRgb5+fm246gQERUVxZgxYxgwYBb16nXF5drPqlUPsWrVozqNWRDTYvRDeXk5Y8aMwRjDiBEjSE1NtR1JhZjk5Kb07TuFzp3HEh1di/37N+oppEFMP2P0w6xZs1ixYgUNGzZkwoQJtuOoECUSRevWV9Ogwdm4XLnExnr2VJeVlXDoUD4JCfUtJ1SH6RbjSeTm5nLvvfcC8Pzzz1O7duA/OFeRJSmp8VEnBXz33UwWLvwDP/zwL4wpt5hMHabFeBIPPvggubm5DBw4kGHDhtmOo8KMMeUUFu6hrKyI9eufZfnyOygs/NF2rIinxXgCO3fu5NVXXyU2NpYXX3zxqJlVlAoEkSh69nyczMwJxMXV4Zdf1rJo0Qhycmbr1qNF+hnjCTRr1owVK1awZs0aOnToYDuOClMiQpMm55OW1oONG6ewa9fnbNw4hd27F9O79zPExibajhhxtBhPIjMzk8zMTNsxVASIj69Djx6P0LjxQNavf5a4uNrExCTYjhWRtBiPIScnhy1btnDppZfajqIiUKNG/ahXryvGuCs+viko2EVUVDSJiY0sp4sM+hmjD2MMo0eP5rLLLmPq1Km246gIFReXSnz8aQCUl5exevUEFi0ayY4d8zHGWE4X/rQYfcyaNYsFCxZQr149rr32WttxlMLtdpGY2JCysiLWrn2KVase1rNmHKbFWMmePXu46667AJg8eTL16+sBt8q+w5Phduv2EDExSezZk8WiRSPZt2+17WhhS4vRyxjDmDFjOHDgAIMHD+b666+3HUmpCiJC8+YXM3Dga9St25mSkl9YvvxONm9+xXa0sKTF6DV79mzef/99kpOTeemll/SYRRWUEhMb0bfvFNq3/xMi0cTH17EdKSzpXmk8W4vPPPMMAM888wzNmjWznEip44uKiqFdu5E0anQuKSmtKu7/+eefadBAZ5QPBN1ixPM2ZeHChbzwwgv8+c9/th1HKb+kpraueGeTk5NDu3btGDduHIcOHbKcLPRpMXqlpKRwxx136FtoFZK++uorCgsLef755xk4cCA//qjnW1eHY8UoIs1EZJGIbBaRTSIy1qmxTtWuXbsYN24cRUVFtqMoVS033HADWVlZNG3alGXLltG9e3cWLVpkO1bIcnKLsQwYZ4zpAPQGxohIRwfHq5Ly8nJGjBjB888/XzGtmFKhrE+fPnzzzTecd9557N27lwsuuIAXXnhBDwg/BY4VozFmjzHmG+/3B4HNQBOnxquqyZMns3DhQtLT0xk/frztOEoFRHp6Op999hn3338/breb++67j++//952rJBTI3ulRaQl0A34+hiP3QzcDNC8efOaiEN2djYPPPAAADNnztQ9eSqsREdH89RTT9GjRw/y8vJo166d7Ughx/FiFJFk4D3gDmPMb64iZYyZAcwAyMzMPKVt/iFDlvr93LKyErKybsLlctGixe+YOfM0Zs488fJ6aUoVit58syHQkHnzPL/fP/+8gpiYROrV6xqQ9YfzvwtHi1FEYvGU4lvGmKC44vjGjX/n4MEfSE5uQadOt9mOo1SNKCz8kdWrJ+B2u+jS5S69LvpJOLlXWoCZwGZjzPNOjVMV5eVllJUVERUVR2bmBJ3rTkWMhIQGtGhxGca4Wb/+GbKzp1BeXmY7VtBycouxL/AHIFtE1nnve9AYM9/BMU8oKiqGHj0eoaBgBykpLWzFUKrGRUXF0KnTbaSktGb9+mfZtm02BQXbycz8n4qrFaojnNwrvdQYI8aYLsaYM71fVkrR7S6lrKwE8JzloqWoIlWLFpfSt+9k4uLqsHfvSrKyRuvFt44hIs582bhxKllZf+bgwR22oyhlXb16XejffwYpKa0oKtqDy3XAdqSgE/aTSOzatYAffnifqKhY3O5i23GUCgqJiY3o1286Bw5soW7dDNtxgk5YbzHm5W1l3bpJAHTqdBt16ujxXEodFhubRHp6j4rbu3cv5ttvX9bLthLGW4ylpfmsXPkQbreLZs0upmXLIbYjKRW0SksPsnbtU5SVFVFYuMs7W3gt27GsCcstRmPcrFnzGEVFu6ld+wy6dr1bZ81R6gTi4lLo2fMJYmKS2L17McuW3U5JyS+2Y1kTlsW4Z88y9u79mri42vTq9STR0fG2IykV9OrX70m/ftNJTGzEgQObycr6M3l5W23HsiIsi7Fx43Pp0uUuevZ8nMTEhrbjKBUyUlNbce65L1O3bieKi/eyZMlf+Pnn30xxEPbCqhgrT6/UqtWVpKV1s5hGqdAUH38aZ5/9N5o2vQAwEXldmbApxsLCH1my5BY9VlGpAIiOjqd79/H07/+/Rx3NESl7rMOiGEtL8/nqq3v59ddNfPfdq7bjKBUWPGeJtay4vXPnZyxbNjYiDggP+WIsKSlh5coHKSjYQWpqG7p2vcd2JKXCjttdypYtr5Cbu46srJvJz99mO5KjQroY3W43f/jDH8jNXU+tWmn07v00sbFJtmMpFXaio+M455wXqVOnPUVFe8jKGs2//vUv27EcE7LFaIzh9ttvZ86cOcTEJNO79zMkJNS3HUupsJWQkM4550ylSZNBuN3FDBkyhMcee4zy8vD73DFki3Hp0qVMmzaN+Ph4zjrrKWrXbms7klJhLzo6nh49HqFjx9GICI8++ijjxo2zHSvgQrYY+/Xrx4wZM3j77bdJSzvTdhylIoaIcPrpN/DRRx/RpEkTRo0aZTtSwIVsMQLcdNNNXHnllbZjKBWRBg8eTE5ODhkZR2bnWbVqlcVEgRPSxaiUsis+/sjptq+99hq9evVi7NixlJaWWkxVfVqMSqmAcLlcxMbGMmXKFPr168e2baF7SI8Wo1IqIEaPHs2SJUto3rw5K1eupFu3brz77ru2Y50SLUalVMCcddZZrFu3jqFDh5Kfn8+wYcMYOXIkBQUFtqNViRajUiqgTjvtNObMmcO0adOoVasWq1atIjo62nasKgnbGbyVUvaICLfccgsDBgzA7XaTkOC5hnteXh6xsbEkJiZaTnhiusWolHJMhw4d6NSpU8XtW2+9la5du7J48WJ7ofygxaiUqhEHDx5k3bp1bN26lYEDB/LHP/6RX34JzssnaDEqpWpESkoKq1ev5rHHHiM+Pp5Zs2bRrl07pk+fjtvtth3vKFqMSqkaEx8fz/jx49mwYQPnnXce+/fv5y9/+Qu9evXC5XLZjldBi1EpVePOOOMMFixYwOzZs2nevDm9evU66iwa27QYlVJWiAhXX301mzdvZtKkSRX3z507l+uuu47Nmzdby6bFqJSyKjExkdq1a1fcnjhxIu+88w4ZGRkMGzaMtWvX1ngmLUalVFB5//33GT16NDExMbz77rt0796diy66iE8++aTGJsXVYlRKBZVmzZoxffp0cnJyuPPOO0lKSuKzzz7jkksu4a233qqRDFqMSqmg1KxZM55//nl27NjBpEmT6Nq1K0OHDq2RsbUYlVJBrW7dutx3332sXbuWpKSaudidniutlAqIIUOWOrr+Dz44x9H1V6ZbjEop5UOLUSmlfGgxKqWUDy1GpZTy4WgxisjFIvKdiGwVkfudHEsppQLFsWIUkWjgReASoCNwnYh0dGo8pZQKFCe3GHsBW40x24wxpcA7wBUOjqeUUgHh5HGMTYCdlW7vAs7yfZKI3Azc7L1ZICLfOZjplIj85q40wNGph48xZiAcN7dD4x1XFcer9utt8edz/HfFZ7xAOWnuIP2dOV7uFlUZy8liPNaPYX5zhzEzgBkO5gg4EVltjMm0naOqNHfNC9XskZ7bybfSu4BmlW43BXY7OJ5SSgWEk8W4CjhdRFqJSBwwDPjQwfGUUiogHHsrbYwpE5FbgU+BaOBVY8wmp8arYSH11r8SzV3zQjV7ROcWY37zsZ9SSkU0PfNFKaV8aDEqpZQPLcZKTnYKo3hM8T6+QUS6V3rsBxHJFpF1IrK6ZpP7lb29iKwQEZeI3F2VZZ1UzdzWXnM/ct/g/R3ZICLLRaSrv8s6rZrZg/k1v8KbeZ2IrBaRc/xd9jeMMfrl+Zw1GsgBWgNxwHqgo89zBgP/xnOMZm/g60qP/QCkBXH2+kBP4Eng7qosG4y5bb7mfuY+GzjN+/0lh39XbL7e1c0eAq95Mkf2m3QBtpzqa65bjEf4cwrjFcDrxuMroI6INKrpoMdw0uzGmL3GmFXAoaou66Dq5LbJn9zLjTG/em9+hec4Xr+WdVh1stvkT+4C421CIIkjJ5RU+TXXYjziWKcwNqnCcwzwmYis8Z7mWJP8ye7EstVV3bFtveZVzf0nPO80TmXZQKtOdgjy11xErhSRLcDHwP+ryrKV6TVfjvDnFMYTPaevMWa3iNQHPheRLcaYrIAmPD6/Tr90YNnqqu7Ytl5zv3OLyEA85XL48y6br3eVxj9Gdgjy19wY8z7wvoicCzwODPJ32cp0i/EIf05hPO5zjDGH/9wLvI9n872mVOf0S5unblZrbIuvuV+5RaQL8ApwhTEmtyrLOqg62YP+NT/MW9ZtRCStqsseXoF+eT6WiAG2Aa048gFths9zLuXonS8rvfcnASmVvl8OXBxM2Ss9dwJH73zxe9kgy23tNffzd6U5sBU4+1R/5iDMHuyveVuO7HzpDvzo/bda5de8Rv4yQuULz17n7/HswXrIe99oYLT3e8Ez+W4OkA1keu9v7X2x1wObDi8bZNkb4vmfMx844P0+9XjLBntu26+5H7lfAX4F1nm/Vp9o2VDIHgKv+X3eXOuAFcA5p/qa6ymBSinlQz9jVEopH1qMSinlQ4tRKaV8aDEqpZQPLUallPKhxRjmRKSgGstefngmEhEZUvm64CLymIgMCkRGnzEXi8gJL2bkz3P8HOsOEUmsdHu+iNSp7np9xhggInkiMv8kz1skIgWB+LlU9WkxquMyxnxojJnkvTkE6FjpsUeMMQusBAucO4CKYjTGDDbGHHBgnCXGmMEneoIxZiBQ49PVqWPTYowQIhIlItNEZJOIfOTdOrra+9gPIvI/IvKNd6699t77R4rIVBE5G7gceMY7110bEZlVafnzRWStd9lXRST+ROv1yZUgIu9459F7F0io9NiF3rkYvxGR2SKSfIzlp3vn3tskIv9TKc/7lZ5zgYjM9VnudqAxsEhEFlXKmyYiLUVki4i8IiIbReQtERkkIstE5D8i0sv7/CTvz7vK+/OfdJYcEWkkIlne13GjiPQ72TKq5mkxRo6hQEugMzAK6OPz+C/GmO7AdOCoCWGNMcvxXOHxHmPMmcaYnMOPiUgtYBbwe2NMZzynX93iz3q9bgGKjDFd8My52MO73jTgYWCQd/nVwF3HWP4h47mOcBegv/cc34VABxFJ9z7nj8BrPj/TFDznyw70bq35agtM9q63PXA9nskU7gYePDw2sNAY0xMYiOc/jqRjrKuy64FPjTFnAl3xnKWhgowWY+Q4B5htjCk3xvwELPJ5/PAW1Ro8BeqvdsB/jTHfe2//Azi3Cus9F3gTwBizAdjgvb83nrfuy0RkHTACaHGM5a8VkW+AtUAGnglIDfAGMNz7mWEfjp46yx//NcZkG2PK8Zxm9oV3vdmVfo4Lgfu9+RYDtfCcZ3wiq4A/isgEoLMx5mAVc6kaoNOORY5jTb1Umcv7p5uq/V4EYr3HOi9VgM+NMdcdd2CRVni24HoaY34VkVl4ygk8W4jzgBI8/yGUnSTn8XIDlFe6Xc6Rn0OAq4wx3/m7UmNMlndKrEuBN0TkGWPM61XMphymW4yRYylwlfezxgbAgCoufxBIOcb9W4CWItLWe/sPwJdVWG8WcAOAiHTC89YVPDNH9z28XhFJFJEzfJZNBQqBPO/PdMnhB4xneqzdeN6Oz6riz+SvT4HbRES8GbudbAERaQHsNcb8LzATzywwKshoMUaO9/DMTLMReBn4GsirwvLvAPd4dzK0OXynMaYEz2d4s0UkG88W1UtVWO90IFlENgD3Aiu9690HjATe9j72FZ7P+ioYY9bjeQu9CXgVWOaz7reAncaYb48z9gzg34d3vpyCx4FYYIOIbPTePpkBwDoRWQtchedzTBVkdHadCCIiycaYAhGph6eA+no/bwxLIjIVWGuMmWkxwwA880he5sdzF3ufq4ftWKafMUaWj7w7I+KAx8O8FNfgeZs9znKUUqCTiMw/0bGM3q3W1gTXRb8ilm4xKqWUD/2MUSmlfGgxKqWUDy1GpZTyocWolFI+tBiVUsrH/wdIh6jNGHujYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[5,4])\n",
    "n, bins, patches = plt.hist(x=IgnTime, bins='auto', color='#0504aa',\n",
    "                            alpha=0.7, rwidth=0.85,density=True)\n",
    "plt.plot(x, gkde.evaluate(x), linestyle='dashed', c='black', lw=2)\n",
    "plt.xlabel('Ignition delay time [s]')\n",
    "plt.ylabel('pdf')\n",
    "# plt.ylim([0,4])\n",
    "plt.savefig('H2_IgniDelayPdf.pdf', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "del(tchem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "pytchem.finalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
